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DETERMINATION  OF  THE  ACOUSTIC  PROPERTIES  OF  FROZEN  SOILS 


I.  ULTRASONIC  VELOCITIES  OF  DILATATIONAL  HAVES  IN  FROZEN  SOILS 

ti 

by 

*  K.  Knuth,  M.  Smith,  R.  Martin  and  Y.  Nakano 


Introduction 

Measuring  the  variation  of  the  acoustic  properties  of  solids  under  variable  physical  conditions 
has  become  well  established  as  an  effective  method  for  investigating  the  physical  structure  of 
solids. 

Although  the  acoustic  properties  of  metals,  plastics,  and  unfrozen  earth  materials  have  been 
widely  explored,  little  attention  has  been  directed  toward  frozen  earth  materials  in  the  past.  Re¬ 
cently  there  has  been  considerable  interest  in  frozen  earth  because  of  military  applications,  such  as 
seismic  monitoring  and  personnel  sensor  detection  in  cold  environments,  and  construction  engineering 
applications,  such  as  the  Trans-Alaskan  pipeline  construction. 

One  of  the  most  important  acoustic  properties  of  solids  is  the  velocity  of  dilatational  waves. 

A  summary  of  dilatational  wave  velocity  data  obtained  in  permafrost  regions  was  compiled  by 
Barnes  (19631.  Eykov  (1966)  reviewed  Russian  works.  Several  laboratory  studies  on  the  subject 
have  been  reported  (Frolov,  1961;  Mullet,  1961;  Kaplar,  1963;  Desai  and  Moore,  1967;  Timur,  1968). 
Frolov  (1961)  measured  velocities  of  30  kHz  dilatational  waves  in  four  different  types  of  frozen 
soils  (sand,  clay,  sandstone  and  silt)  in  the  temperature  range  from  -20°C  to  20°C.  Muller  (1961) 
measured  the  velocity  in  water-saturated  sand  and  clay  of  various  porosities  as  a  function  of 
decreasing  temperature.  His  results  indicate  that  with  increasing  ice  content  the  velocity  decreases 
for  sand  and  increases  for  clay.  Kaplar  (1963)  measured  both  dilatational  and  shear  wave  velocities 
in  various  frozen  earth  materials  in  the  temperature  range  from  0°C  to  -20°C  by  the  resonant  bar 
method,  in  which  either  flexural,  longitudinal  or  torsional  vibrations  were  induced  by  electromagnetic 
means. 

Recently  Timur  (1968)  measured  dilatational  wave  velocities  in  various  earth  materials  between 
26°C  and  -36°C  by  the  pulse  first-arrival  technique,  in  which  the  time  required  for  an  elastic  wave 
to  traverse  a  sample  of  known  length  is  determined.  He  measured  velocities  with  both  descending 
and  ascending  temperature  and  found  that  the  two  measurements  generally  do  not  agree,  the  degree 
of  discrepancy  depending  on  the  specimen.  It  has  recently  been  shown  that  freezing  and  thawing 
bring  about  a  dramatic  redistribution  of  water  and  a  reorientation  of  particles,  particularly  in  fine¬ 
grained  earth  media  such  as  clay  (Anderson  and  Hoekstra,  1965a  and  1965b;  Anderson  and  Tice, 

1970).  It  is  possible  to  consider  a  hysteresis  of  velocity  reported  by  Timur  (1968)  as  a  result  of 
such  structural  change  in  the  specimen. 

This  report  covers  the  first  phase  of  an  investigation  of  the  relationship  between  acoustic 
properties  of  frozen  soils  and  soil  structure  as  well  as  constituents.  The  velocities  of  dilatational 
waves  in  three  standard  soils  were  measured  with  the  pulse  first-arrival  technique.  A  hysteresis  of 
velocities  similar  to  that  obtained  by  Timur  was  observed. 
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Experimental  Procedure 

Sample  preparation 

■  Three  standard  types  of  soils.  20-30  Ottawa  sand,  Hanover  silt  and  Goodrich  clay,  were 
tested  under  fully  water-saturated  conditions.  Figure  1  shows  the  gradation  curves  obtained 
according  to  ASTM  (American  Society  for  Testing  and  Materials)  test  ptocedures.  We  prepared 
circular  cylindrical  samples,  either  2.5  cm  diameter  x  about  15  cm  length  oi  2.5  cm  diameter  x 
about  90  cm  length.  To  prepare  sand  or  silt,  dry  soil  was  first  packed  into  1-inch  Tygon  tubing 
encased  in  a  copper  jacket  and  was  tamped  or  vibrated  until  a  specified  dry  density  was  attained. 
Then  water  was  sucked  into  the  sample  by  the  use  of  vacuum.  To  prepare  clay,  water-saturated 
clay  was  packed  into  Tvgoii  tubing  in  orde"  to  maintain  uniform  density  throughout  the  sample. 
After  the  ends  of  the  tubing  were  sealed  with  aluminum  pings,  the  sample  protected  by  the  tubing 
and  the  copper  jacket  was  frozen.  When  the  sample  was  ready  for  testing,  the  copper  jacket  was 
removed  and  the  plugs  were  replaced  with  transducer  assemblies. 


Figure  I.  Gradation  curves  (or  20-30  Ottawa  sand,  Hanover  silt,  and  Good¬ 
rich  clay.  Gs  specific  gravity. 


Temperature  control 

Copper-constantan  thermocouples  were  inserted  to  the  center  of  the  sample  to  monitor  tempera¬ 
ture.  To  maintain  a  constant  temperature  during  the  experiment,  a  Forma  Scientific  R  idel  2095 
bath  was  used.  The  sample  was  placed  in  a  cooling  jacket,  through  which  the  cooling  fluid  was 
circulated  by  a  puiup  via  the  bath.  The  temperature  of  the  sample  was  kept  constant  within  ±0.1  °C. 

Velocity  measurement 

The  velocities  of  propagating  waves  in  the  sample  were  measured  with  the  pulse  first-arrival 
technique  (Kolskv,  1963).  At  either  end  of  the  sample  a  transducer  of  0.5-inch-diameter  x 
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0.25-inch-thick  PZT4,  which  was  bonded  to  an  aluminum  disk  1  inch  in  diameter  and  0.25  inch 
thick  with  silver  epoxy,  was  attached  to  the  surface  of  the  sample  with  a  few  drops  of  silicone 
oil  for  better  coupling.  One  of  the  transducers  served  as  a  transmitter  and  the  other  as  a  receiver. 

A  Hewlett  Packard  Model  214A  pulse  generator  supplied  pulses  of  about  1  second  duration  to 
the  transmitter  with  a  repetition  rate  of  500  to  1000  pulses  per  second.  The  receiver,  which  was 
connected  to  a  Tektronix  50A  dual-trace  oscilloscope  via  a  Krohn-Hite  Model  3202  filter,  displayed 
the  received  signal.  A  filter  was  used  in  a  bandpass  mode  passing  10  kHz  to  1  MHz.  The  other 
oscilloscope  trace  was  used  for  a  Computer  Measurements  Company  Model  All  digital  time  delay 
generator  providing  accurate  measurements  of  arrival  time  with  an  error  of  less  than  +0.1  sec.  Both 
traces  on  the  oscilloscope  were  triggered  by  the  pulse  generator.  This  system  was  checked  by  the 
use  of  a  standard  medium,  such  as  water,  copper  or  aluminum,  prior  to  measurements  on  frozen  soils. 
All  tests  gave  velocities  within  1%  of  handbook  values. 

Accuracy  of  the  pulse  method 

The  behavior  of  elastic  waves  in  any  bounded  medium  necessarily  entails  various  effects  due 
to  the  presence  of  the  boundaries.  The  effects  of  boundaries  in  circular  cylinders  have  been  well 
studied  (Kolsky,  1963).  Figure  2  depicts  the  theoretical  group  velocity  of  the  first  six  axially 
symmetric  modes  in  an  aluminum  rod  1.5  inches  in  diameter.  The  calculations  were  performed 


Figure  2.  Group  velocity  versus  frequency  for  a  circular  cylinder  1.5  in.  in  diameter,  having  a  com¬ 
pressions!  velocity  of  6.42  km/sec.  The  expected  “bar”  velocity  is  3.04  km/sec.  The  crosses 
are  observed  group  velocities  in  a  1.5  in.  aluminum  rod,  327.5  cm  long. 
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using  handbook  values  for  shear  and  dilatational  velocities.  The  lowest  [node  approached  the  bar 
velocity  in  the  low-frequency  limit.  It  is  clear  from  Figure  2  that  the  mode  structure  of  a  cylinder 
is  fairly  complex.  The  crosses  represent  observations  of  group  arrival  versus  frequency  taken, 
with  our  apparatus,  on  a  327.5-em  aluminum  rod,  1.5  inches  in  diameter.  The  primary  uncertainty 
in  the  data  lies  in  the  assigned  frequency,  which  was  obtained  by  a  period  measurement  of 
adjacent  peaks.  The  agreement  is  quite  satisfactory  within  the  limits  .imposed  by  this  error.  The 
spurious  points  at  about  80  kHz  may  be  attributable  to  mode  coupling  at  the  rod's  support  points 
although  we  did  not  attempt  to  verify  this. 

As  is  seen  in  Figure  2,  successively  higher  modes  possess  group  velocity  maxima  which 
increase  in  both  frequency  and  velocity.  In  a  general  way,  these  maxima  approach  the  dilatational 
velocity  of  an  elastic  medium  increasingly  closely.  For  example,  the  14th  mode  possesses  a  maxi¬ 
mum  at  1.15  MHz  and  a  velocity  of  97%  of  dilatational  velocity. 

In  view  of  these  calculations  dilatational  wave  travel-time  measurements  should  be  made  at  the 
highest  feasible  frequency.  Care  should  be  taken  to  ensure  operation  in  a  region. where  the  group 
velocity  is  close  enough  to  the  dilatational  velocity  to  achieve  the  desired  accuracy. 

So  far  we  have  discussed  the  accuracy  of  the  method  applied  for  elastic  solids.  Since  any 
real  material  deviates  from  an  ideal  elastic  solid  in  one  way  or  another,  the  accuracy  of  the  method 
also  depends  on  the  anelastic  behavior  of  the  solids  examined. 


Results  and  Discussion 

The  results  of  the  experiment  are  presented  in  Figures  3-5  where  dilatational  velocities  are 
plotted  as  a  function  of  temperature.  Originally  we  intended  to  evolve  a  technique  that  would 
allow  simultaneous  measurement  of  both  dilatational  and  shear  velocities  using  long  samples  based 


Figure  3.  Dilatational  velocity  vs  temperature  for  Ottawa  sand.  The  differences  in  density 
and  velocity  reflect  the  difference  in  porosity. 
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upon  the  theory  of  guided  elastic  waves  in  a  cylinder.  This  effort  has  been  unsuccessful  due  to 
unexpected  high  attenuation  in  the  samples  tested  For  Ottawa  sand  and  Hanover  silt,  90-cm-long 
samples  barely  allowed  measurement  of  dilatational  velocities  in  the  frozen  state,  and  did  not 
allow  measurement  in  the  unfrozen  state.  Ninety-cm-long  Goodrich  clay  did  not  yield  any  reliable 
measurement  even  in  the  frozen  state.  Despite  the  limited  experimental  data,  we  are  able  to 
present  the  following  observations. 

Velocity  versus  temperature  curves 

The  decrease  in  dilatational  wave  velocity  as  an  initially  frozen  water-saturated  soil  is 
thawed  appears  to  be  a  direct  consequence  of  the  change  in  state  of  the  water.  Consider  the  soil 
as  a  two-component  mixture,  that  is,  a  granular  framework  whose  interstitial  spaces  are  filled 
with  water.  A  high-frequency  wave  traveling  along  any  path  through  the  sample  will  travel  part  of 
the  tune  through  the  crystalline  framework  and  part  of  the  time  through  either  water  or  ice,  depend¬ 
ing  on  the  temperature.  Since  the  dilatational  wave  velocity  is  about  3.0  km/sec  in  polycrystalline 
ice  but  only  1.49  km/sec  in  water,  we  expect  that  the  travel  time  along  the  same  path  in  a  frozen 
sample  will  be  less  than  the  corresponding  travel  time  in  the  thawed  sample.  The  velocity  of  a 
wave  propagating  in  the  crystalline  framework  is  essentially  constant  over  this  temperature  range. 
Hence,  the  observed  velocity  for  frozen  soil  should  be  greater  th\n  for  the  same  sample  thawed. 

In  light  of  this  explanation  of  the  variations  of  dilatational  velocity  in  saturated  soil  as  a 
function  of  temperature,  it  might  be  asked  whether  or  not  the  obser/ed  change  in  velocity  is  in 
some  way  proportional  to  the  amount  of  water  in  the  sample.  This  suggests  that  perhaps  each 
component  contributes  to  the  observed  slowness  in  proportion  to  its  relative  abundance  in  the 
sample  and  the  average  compress  ional  wave  velocity  of  the  individual  components.  Averaging 
techniques  for  two-component  systems  have  been  employed  successfully  to  obtain  compressibilities 
and  moduli  of  many  minerals  which  are  available  only  in  small  quantities  or  in  finely  divided 
particles  not  suitable  for  bulk  testing.  They  are  first  mixed  with  an  isotropic  material  with  known 
elastic  constants.  The  unknown  elastic  parameter  is  then  computed  from  the  average  properties 
of  the  composite  material  (Anderson,  1963;  Chung  and  Buessman,  1967;  Brace  et  al.  1969).  The 
most  notable  method  is  to  use  Voigt  and  Reuss  averages  to  obtain  upper  and  lower  bounds  and  then 
select  a  composition  for  which  the  spread  of  the  bounds  is  a  minimum.  These  methods  however 
are  not  directly  applicable  in  a  straightforward  way  to  the  simple  averaging  of  dilatational  wave 
velocities  for  a  two-component  system. 

Timur  (19C9)  measured  the  dilatational  wave  velocity  in  water-saturated  porous  sandstone  as 
a  function  of  temperature.  He  observed  a  change  in  velocity  similar  to  the  one  we  observed  as  the 
temperature  of  tbo  sample  was  increased  from  -24°C  to  +24°C.  Moreover,  only  an  insignificantly 
small  decrease  in  velocity  was  observed  on  a  dry  sandstone  sample  as  the  temperature  was  raised 
from  Delow  the  fpeezing  point  of  water  to  room  temperature.  These  observations  also  strongly 
suggest  that  the  change  in  dilatational  velocity  as  a  frozen  sample  is  thawed  is  attributable  solely 
to  the  state  of  water  in  the  pore  spaces  of  the  rock.  Using  an  argument  similar  to  the  one  we  out¬ 
lined  above  to  explain  the  variation  in  velocity  as  our  samples  thawed,  Timur  (1968)  proposed  a 
simple  time -averaging  method  to  compute  the  theoretical  velocity  of  a  two-component  system  based 
on  the  percentage  of  each  component  in  the  sample  and  the  velocity  of  that  component.  This 
technique  assumes  that  the  tiavel  time  for  a  dilatational  wave  through  the  sample  is  the  travel 
time  for  each  component  computed  according  to 
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where  vQ  is  the  observed  velocity  for  the  sample,  Vj  and  v2  are  the  velocities  of  each  component, 
and  A  is  the  relative  volume  of  one  component.  For  a  water-saturated  sample,  the  relative  volume 
of  water  present  is  equal  to  the  porosity  of  the  sample. 

Applying  eq  1  to  his  results,  Timur  found  good  agreement  with  his  experimental  results. 
Typically  the  observed  velocity  fell  within  5%  of  the  predicted  value.  Our  results  on  20-30  Ottawa 
sand  with  a  porosity  of  37.5%  (Fig.  3)  were  compared  with  the  predicted  velocities  obtained  using 
the  time-averaging  equation.  For  the  frozen  soil  the  predicted  and  observed  values  agreed  to 
better  than  2%.  However,  for  saturated  soil  above  0°C  agreement  was  poor.  The  observed  velocity 
was  2.78  km/sec  whereas  the  predicted  value  using  time-averaging  was  not  immediately  applicable 
4  to  the  silt  sample  because  the  composition  and  relative  abundance  of  the  solid  particles  were  not 

known. 

The  inapplicability  of  the  time-averaging  technique  to  unfrozen,  saturated  soil  presents  a 
difficulty.  Why  does  the  time-averaging  method  appear  to  be  satisfactory  for  both  frozen  and  un¬ 
frozen  porous  rocks  but  only  to  frozen  Ottawa  sand?  The  fact  that  the  compressional  wave  velocity 
is  less  than  the  anticipated  value  based  on  time-averaging  in  thawed  soil  but  in  good  agreement  for 
frozen  soil  suggests  that  the  compressional  velocity  is  largely  determined  by  the  compressibility 
of  the  interstitial  water  rather  than  by  the  compressibility  of  the  mineral  solids.  It  appears  then 
that  for  a  consolidated  and  lithified  elastic  continuum,  such  as  a  porous  rock,  the  solid  and  liquid 
can  be  time-averaged  in  direct  proportion  tc  their  relative  volumes  to  obtain  a  realistic  value.  In 
the  case  of  a  non-lithified,  unfrozen  soil,  however,  time  averaging  does  not  seem  to  apply  due  to 
the  discontinuous  nature  of  the  mineral  grains. 

Hamilton  (1970)  studied  the  velocity  of  water-saturated  marine  sediments  as  a  function  of 
porosity  and  grain  diameter.  He  found  that  the  velocity  typically  increased  from  about  1.50  km/sec 
to  1.86  km/sec  as  the  porosity  was  decreased  from  80  to  30%.  A  similar  increase  in  velocity  was 
observed  when  the  grain  diameter  was  increased  from  1  tc  1000/i.  Applying  the  time-averaging 
equation  to  Hamilton’s  results  for  sand  produces  no  agreement  between  the  predicted  and  observed 
results.  Our  results  and  those  obtained  by  Hamilton  suggest  that  time  averaging  is  not  an  applicable 
method  for  predicting  velocities  in  water-saturated  soils.  The  breakdown  of  the  time-averaging 
approach  to  velocities  when  the  solid  minerals  do  not  form  an  interconnected  framework  suggests 
that  further  justification,  on  a  sound  physical  basis,  is  required  before  time  averaging  can  be 
accepted  wholeheartedly  even  for  frozen  soils. 

Ideally  we  would  like  to  have  a  single  theory  that  would  predict  the  velocity  in  both  frozen  and 
unfrozen  soils  using  easily  measured  properties  of  the  components  such  as  relative  abundance, 
mineral  composition,  velocity  of  the  minerals,  and  other  readily  obtainable  elastic  properties. 

»  Presently  no  such  theory  exists.  As  a  major  protion  of  our  future  research  effort  on  the  geophysical 

properties  of  frozen  soils,  this  problem  will  be  analyzed  in  detail  both  empirically  and  analytically 
to  achieve  a  workable  relation  between  porosity,  mineral  composition,  water  content  and  sample 
•  velocity. 

Hystei  esis  In  the  velocity  during  a  freeze-thaw  cycle 

It  is  evident  that  a  strong  correlation  exists  between  dilatational  velocities  and  unfrozen  water 
content.  Then  it  might  be  asked  whether  or  not  the  observed  hysteresis  in  the  velocity  during  a 
freeze-thaw  cycle  is  also  caused  by  the  hysteresis  of  unfrozen  water  content.  The  low- temperature 
phase  composition  of  interfacial  water  has  been  a  topic  of  continuing  interest.  Ncrsesova  and 
Tsytovich  (1963)  conducted  experimental  investigations  to  determine  the  phase  composition  of  water 
in  various  frozen  soils  by  calorimetric  methods.  Unfrozen  water  contents  obtained  by  them  for 
typical  non-saline  soils  are  shown  in  Figure  6.  In  granular  soil,  pores  are  comparatively  large  and 
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Figure  6.  Unfrozen  water  contents  in  typical  non¬ 
saline  soils:  1)  quartz  sand,  2)  sandy  loam,  3)  loam, 

4)  clay,  and  5)  clay  containing  montmorillonite. 

almost  all  of  the  water  freezes  at  the  freezing  point  of  water.  However,  clay  and  silt  have  fine 
pores,  in  which  a  significant  portion  of  the  water  remains  unfrozen  in  a  liquid  or  semiliquid  state. 
Those  observations  are  consistent  with  the  velocity-temperature  relation. 

Several  studies  concerning  phase  composition  have  been  reported.  However,  the  hysteresis  of 
phase  composition  during  a  freeze-thaw  cycle  has  not  been  discussed  explicitly,  because  equilibrium 
phase  composition  has  been  the  main  subject.  Anderson  and  Hoekstra  (1965a,  1965b)  have  recently 
shown  that  freezing  and  thawing  bring  about  a  dramatic  redistribution  of  water  and  a  reorientation 
of  particles  in  clay.  They  studied  the  changes  in  apparent  d  (001)  spacing  in  Wyoming  bentonite 
dining  the  freeze-thaw  cycle  by  X-ray  diffraction.  They  found  the  hysteresis  loop  in  the  spacing 
during  the  freeze-thaw  cycle  due  mainly  to  supercooling.  However,  if,  during  cooling,  the  sample 
was  nucleated  artificially  the  d  (001)  spacing  dropped  immediately  and  the  cooling  and  warming 
curves  nearly  coincided.  The  d  (001)  spacing  indicates  the  amount  of  interlamellar  water  consti¬ 
tuting  the  gel  structure.  Freezing  altered  the  gel  structure  and  most  of  the  interlamellar  water 
was  expelled  on  complete  freezing;  as  a  result,  ice  must  form  in  extralamellar  regions.  If  the 
interlamellar  water  corresponds  to  the  ’‘unfrozen"  water  (Anderson  and  Hoekstra,  1965a,  1965b), 
the  hysteresis  of  unfrozen  water  can  occur  most  probably  in  laboratory  experiments,  where  no  control 
on  nucleation  is  made.  Since  almost  no  interlamellar  water  exists  in  granular  sand,  such  as  Ottawa 
sand,  it  is  consistent  that  we  did  not  observe  anv  hysteresis  in  velocities  for  Ottawa  sand. 
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Although  the  unfrozen  water  content  is  one  of  the  most  important  variables  in  determining 
dilatational  velocities  in  fine-grained  frozen  soils,  determining  the  amount  of  unfrozen  water  is 
quite  elaborate.  Efforts  have  been  initiated  to  evolve  a  technique  allowing  simultaneous  deter¬ 
mination  of  acoustic  velocities  and  unfrozen  water  contents. 
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n.  DETERMINATION  OF  A  LINEAR  VISCOELASTIC  CONSTITUTIVE  EQUATION  FOR 
FROZEN  EARTH  MATERIALS  BY  THE  USE  OF  THE  RESONANCE  COLUMN  TECHNIQUE 

by 

Y.  Nakano  and  H.  Stevens 


Introdnctloo 

During  the  past  decade  considerable  emphasis  in  soil  mechanics  research  has  been  placed  on 
the  behavior  of  soils  under  dynamic  loading.  This  behavior  depends  strongly  on  the  nature  of 
loading:  stress,  stress  rate,  frequency,  etc.  Consequently,  various  experimental  techniques  have 
been  developed  to  determine  the  dynamic  behavior  of  soils  subjected  to  a  specific  kind  of  loading. 
Among  these,  the  resonance  column  technique  has  been  extensively  used  to  investigate  unfrozen 
soils  under  relatively  weak  harmonic  loading  (stress  less  than  100  psi)  in  the  frequency  range  of 
10'1  -  104  Hz  (Hardin  and  Richart,  1963;  Hardin  and  Mossbarger,  1966;  Hardin  and  Black,  1968; 
Hardin  and  Drnevich,  1970). 

In  the  resonance  column  method  a  cylindrical  column  of  material  is  subjected  to  a  steady 
sinusoidal  loading,  rather  in  the  torsional  or  longitudinal  mode.  When  a  specimen  of  soil  is  under 
vibrational  loading,  the  stress-strain  relation  creates  a  hysteresis  loop.  Two  parameters  have  been 
used  to  dtfine  this  relation  (Hardin  and  Dmevich,  1970).  These  parameters  are  the  modulus  defined 
by  the  slope  of  a  line  through  the  ends  of  the  loop,  and  the  area  of  the  loop,  which  is  a  measure  of 
the  damping  capacity  of  the  soil.  Another  way  of  defining  the  stress-strain  relation  is  to  determine 
the  complex  modulus  according  to  linear  viscoelastic  theory  (Lee,  1963;  Hardin,  1965). 

The  dynamic  behavior  of  frozen  soils  is  less  complex  than  that  of  unfrozen  soils;  but  no 
comprehensive  description  of  either  material  has  been  obtained.  One  useful  and  practical  approach 
towards  a  quantitative  description  of  dynamic  behavior  is  determination  of  a  constitutive  relation 
based  upon  mechanics  of  continua.  Stevens  (1967)  used  a  linear  viscoelastic  model  for  interpreting 
resonance  column  experiments  and  determined  complex  shear  and  Young’s  moduli.  In  his  experi¬ 
ments  the  maximum  dynamic  stress  in  the  specimen  was  varied  from  0.1  psi  to  about  5.0  psi,  where 
the  complex  modulus  was  found  to  be  weakly  dependent  on  the  stress  level.  Despite  such  non-linear 
behavior  in  frozen  soils,  the  linear  viscoelastic  constitutive  equation  is  considered  a  good  first 
approximation  under  low  stress  loading.  In  the  present  work,  efforts  were  made  to  determine  the 
simplest  linear  viscoelastic  constitutive  equation  that  can  describe  the  dynamic  behavior  of  frozen 
soils  under  both  torsional  and  longitudinal  vibrations  consistently.  In  practical  applications  we 
encounter  various  types  of  disturbance:  steady  and  unsteady  disturbances,  plane  and  surface  waves, 
and  cylindrical  and  spherical  waves.  Once  the  constitutive  equation  is  determined,  it  is  possible 
to  predict  the  response  of  frozen  soils  subjected  to  such  disturbances. 
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Expt  1  mental  Determinations  of  Coaplex  Modulus 

In  a  linear  viscoelastic  material,  when  the  stress  a  varies  sinusoidally  with  time  at  an  angulat 
frequency  a>  the  strain  <  varies  with  time  at  the  same  frequency  but  there  is  a  phase  lag  8  between 
stress  und  strain.  The  stress  and  strain  relation  for  linear  viscoelastic  solids  is  generally  express¬ 
ed  by  the  following  equation. 

a  -  Et  (1) 

where  E  is  defined  as  the  complex  Young's  modulus  and  is  a  function  of  oj.  The  complex  shear 
modulus  C  is  defined  similarly.  We  describe  a  method  of  determining  these  complex  moduli  in  the 
following. 

Method  of  test 

A  vertical  cylinder  of  soil  is  subjected  to  steady-state  sinusoidal  vibration  in  the  torsional  or 
longitudinal  mode  at  the  lower  base  end  with  the  other  end  free  except  for  a  light,  relatively  rigid 
cap.  The  input  and  output  stress  waves  are  observed  and  measured  by  piezoelectric  accelerometers 
attached  to  the  base  plate  and  cap  plate  at  each  end  of  the  specimen.  The  peak  acceleration  and 
the  frequency  are  recorded.  The  drive  frequency  may  be  any  value  above  the  so-called  “rigid  body 
frequency"  and  within  the  limits  of  the  drive  motors,  if  the  phase  angle  between  input  and  output 
waves  can  be  accurately  measured;  otherwise,  the  specimen  must  be  excited  at  a  known  resonance. 
The  ratio  of  output  to  input  amplitudes  and  the  frequency,  together  with  the  specimen  properties  of 
density  and  length,  are  required  to  compute  the  desired  parameters. 

Apparatus 

The  complete  test  apparatus  includes  a  device  for  holding  the  specimen,  drive  motors  and 
transducers  for  measuring  the  response,  control  and  readout  instrumentation,  and  auxiliary  molds  and 
equipment  for  specimen  preparation.  This  apparatus  has  been  discussed  in  detail  by  Hardin  (1966) 
and  Stevens  (1967).  At  present  the  apparatus  does  not  include  a  pressure  cell,  and  the  specimens 
are  tested  unconfined. 


Computation  of  complex  moduli  and  results 

First  we  consider  longitudinal  vibration.  If  the  wavelength  of  the  standing  wave  is  long  in 
comparison  with  the  diameter  of  the  specimen,  a  rod  condition  is  approximately  true.  The  equation 
of  motion  for  longitudinal  vibration  is; 


„  02u  d2u 
E =  p 

<?x2 

where 

E  complex  Young's  modulus 
p  density  of  the  medium 
t  time 

u  displacement  along  the  coordinate  x 
x  coordinate  (Lagrangian). 


(2) 
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At  the  driven  end  x  a  0  the  system  is  given  a  sinusoidal  displacement,  namely, 

u(0,  t)  =.  uQeiojt  (3) 

At  the  other  end  x  -  L  the  effect  of  the  mass  of  the  cap  resting  on  top  of  the  specimen  is  con¬ 
sidered: 


SE  l)  -  m  d*u(L'  0 
dx  ftZ 


(4) 


wliere 

S  a  cross-sectional  area  of  the  specimen 
m  =■  end  mass. 

Solving  eq  2  with  the  boundary  conditions,  eq  3  and  4,  we  obtain  the  following  relationship  for  the 
ratio  of  bar  end  displacements  (or  accelerations),  z: 


z  ,  =  I  860  PL  1 

lu  (0,  t)  I  II  -  y  tan  pL  I 

where 

£ 

mo)2 

y  ~pSE  * 

In  more  detail  (Norris  and  Young,  1970): 


Re(z  *)  =  cosh  \  (cos£  -  QL£sin£) 

+  Ql  £  tan  ^  cos  £  sin  ^  tan  ^  j 


Im  (z~l) 


sinh  ^  tan  ^  j  (sin  £  +  QL  £  cos  £)  + 


+  QL£tan  |  sin£ cos  ^ftan 


(5) 


(6a) 


(6b) 


where 


8  imlEL 

Re(E) 


14 


DETERMINATION  OF  THE  ACOUSTIC  PROPERTIES  OF  FROZEN  SOILS 


V 


lp 


phase  velocity  = 


8 

sec  - 
2 


<?L 


m 

pSL 


The  condition  of  specimen  resonance  occurs  when  the  ratio  z  is  a  maximum.  As  the  frequency 
is  increased  from  zero,  the  first  maximum  is  the  fundamental  resonance  and  successive  maximums 
indicate  the  harmonics.  For  the  condition  of  resonance,  where  z  is  a  maximum  and  z~8  is  a  mini¬ 
mum, 


(7) 


Equations  5  and  7  may  be  solved  simultaneously  with  a  computer  using  an  iterative  process  for 
the  expressions  and  tan  (S/2).  At  each  resonance  point  we  can  determine  the  complex  Young’s 
modulus  from  observed  values  of  (o  and  z. 

The  equation  for  torsional  vibration  is  perfectly  analogous.  In  the  torsional  mode  QL  is  re¬ 
placed  by  Qt,  which  is  defined  as 


r4p  L 

C^C  c 

r4pL 


(8) 


where 

r  =  radius  of  the  specimen 
rc  radius  of  the  cap 
pc  =  density  of  the  cap 
Lc  =  length  of  the  cap. 

The  dynamic  stress  in  the  sample  varies  along  its  length  as  a  damped  sine  wave,  with  the 
maximum  at  the  node  nearest  the  bottom  or  input  end.  At  resonance,  this  node  is  very  close  to  the 
bottom  plane  of  the  sample  and  the  stress  computed  for  correlative  purposes  is  computed  for  the 
bottom  plane. 

It  is  difficult  or  impossible  to  accurately  control  the  stress  in  the  sample  during  the  test  be¬ 
cause  of  the  resonance  phenomenon.  Closest  control  car.  be  obtained  by  keeping  the  input  accelera¬ 
tion  level  g  constant.  As  g  is  directly  proportional  to  frequency  and  a  wide  frequency  range  is 
required,  it  is  not  very  practical  to  control  g  at  a  relatively  high  value. 

Consequently,  no  attempt  is  made  to  keep  either  stress  or  g  constant.  Instead,  the  drive  force 
is  held  constant  while  the  frequency  is  varied  over  its  entire  range.  Then  the  drive  force  is  in¬ 
creased  an  arbitrary  amount,  and  so  on.  Thus  we  ensure  measurements  over  a  range  of  stress  levels 
without  knowing  in  advance  what  the  values  will  be. 

The  frequencies  at  which  measurements  are  taken  cannot  be  predetermined,  as  only  the 
resonant  condition  is  used  and  resonant  frequency  varies  with  sample  mass  and  stiffness.  The  first 
four  or  five  resonances  are  usually  used. 

It  is  desirable  to  determine  the  moduli  and  loss  angles  for  a  given  frequency  and  a  given  stress 
or  strain,  not  only  because  these  relationships  are  required,  but  to  allow  comparisons  between  values 
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for  different  samples.  To  accomplish  this,  modulus  and  tan  S  are  each  plotted  versus  the  computed 
dynamic  stress  at  one  resonance.  A  smooth  curve  is  drawn  through  the  points  and  the  modulus  of 
tan  8  value  picked  off  for  a  given  stress.  The  resonance  frequency  is  obtained  by  interpolating 
between  the  two  adjacent  measured  frequencies.  If,  say,  four  resonances  are  used,  four  values  at 
four  frequencies  at  a  constant  stress  are  obtained.  A  plot  of  modulus  or  tan  8  versus  frequency  is 
then  prepared  and  a  smooth  curve  drawn  through  these  points.  A  value  for  modulus  or  tan  8  for  a 
given  frequency  and  a  given  stress  can  then  be  obtained.  The  same  process  is  used  to  obtain 
values  for  about  three  stress  levels  and  at  least  three  frequencies. 

Usually,  the  test  measurements  plot  in  such  a  way  that  there  is  little  question  as  to  the  shape 
of  the  curve  drawn  through  the  points.  However,  scatter  does  occur,  particularly  for  the  tan  8 
measurements.  Some  guidelines  are  used  in  drawing  the  curves.  It  is  assumed  that  the  modulus 
decreases  or  is  constant  with  increasing  stress  and  increases  or  is  constant  with  increasing  fre¬ 
quency.  Tan  8  is  assumed  to  increase  with  increasing  stress  but  there  are  no  good  guidelines  for 
the  relationship  to  frequency.  Usually,  the  trend  is  to  a  decrease  with  increasing  frequency  but  at 
times  there  is  a  strong  trend  to  a  maximum  peak  at  a  particular  frequency  within  the  test  range. 

The  results  of  the  experiment  are  presented  in  Table  I.  The  specimens  tested  included  several 
standard  frozen  soils  as  well  as  polycrystalline  ice.  The  following  variables  are  also  listed  in  the 
table: 

L  =  length  of  specimen  (cm) 

D  =  diameter  of  specimen  (cm) 
pw  =  wet  (total)  density  (g/cm1) 

W  «-•  water  content  (g  water/g  dry  soil)  (%) 

PD  -  dry  density  (g/cm1) 

P0  =  porosity  (void  volume/total  volume)  (%) 

VA  -  void  ratio  (void  volume/dry  soil  volume)  (-) 

Sw  =  saturation  (water)  (%) 

Sj  =  saturation  (ice)  (%) 

T  =  test  temperature  (°C) 
t  =  frequency  (kHz) 

E *  =  absolute  value  of  complex  Young's  modulus  (Kbar) 
y]p  =  phase  velocity  of  longitudinal  wave  (km/sec) 

G*  =  absolute  value  of  complex  shear  modulus  (Kbar) 

Ftp  =  phase  velocity  of  torsional  wave  (km/sec) 

=  dynamic  stress  imposed  on  the  bottom  circumference  of  the  specimen  (psi). 

The  values  of  E *  and  G*  in  the  table  are  either  interpolated  or  extrapolated  from  those  at 
resonance  frequencies  in  order  to  show  the  properties  of  different  frozen  soils  at  the  same  frequencies. 
The  phase  velocities,  F,p  and  V(p,  are  computed  based  upon  complex  moduli  E  and  G  respectively. 
Gradation  curves  of  the  soils  obtained  using  ASTM  (American  Society  for  Testing  and  Materials)  test 
procedures  are  also  presented  in  Figure  7. 

Despite  the  limited  data  it  is  possible  to  describe  some  general  trends  in  the  dynamic  behavior 
of  frozen  soils.  Stevens  (1967)  found  several  important  parameters  affecting  such  behavior:  soil 
type,  ice  content,  void  ratio  and  frequency. 
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Table  I  (Coal'd). 
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U.S.  Std.  Sieve  No. 

4  10  40  200 


Coarse 


Curve 

Description 

Specific 

gravity 

Liquid 

limit 

Plastic 

limit 

Plastic 

index 

Non¬ 

plastic 

A 

20-30  Ottawa  sand 

2.65 

X 

B 

-100-200  Ottawa  sand 

2.65 

X 

C 

100  Lebanon  till 

2.86 

X 

D 

Hanover  silt 

2.74 

X 

E 

Fairbanks  silt 

2.70 

X 

F 

Manchester  silt 

2.73 

X 

O 

Suffteld  clay 

2.69 

45.0 

24.0 

21.0 

Figure  7.  Gradation  curves. 


For  ice-saturated  or  almost  saturated  non-plastic  frozen  soils  a  strong  correlation  was  found 
between  modulus  and  void  ratio  (Fig.  8,  9).  Since  such  frozen  soils  have  only  two  constituents, 
soil  minerals  and  ice,  the  moduli  are  bounded  below  by  those  of  ice  and  above  by  those  of  rock. 

The  moduli  of  several  standard  rocks  measured  by  Simmons  and  Brace  (1965)  are  plotted  in  these 
figures  to  show  the  upper  bound. 

Frozen  soil  specimens  have  a  fundamental  resonance  of  the  order  of  2.0  kHz  longitudinally  and 
1.0  kHz  torsionally.  In  this  experiment  soil  was  usually  tested  up  to  the  third  and  fourth  resonance. 
Therefore,  the  range  of  frequency  is  about  1.0  kHz  to  10  kHz.  In  this  range  the  modulus  increases 
and  tan  8  decreases  with  increasing  frequency.  Tan  8  is  approximately  a  reciprocal  of  the  quality 
factor  Q,  which  is  defined  as  the  ratio  of  the  energy  carried  by  a  wave  to  the  energy  dissipated  per 
radius  of  phase  shift  (Thurston,  1964)  and  is  commonly  used  by  seismologists  (Jackson  and  Anderson, 
1970).  For  ice-saturated  frozen  soils  the  correlation  between  Q  and  void  ratio  is  not  so  strong  as 
that  between  modulus  and  void  ratio. 
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Figure  8.  E*  vs  void  ratio,  t  Granite,  Stone  Moun¬ 
tain,  Georgia.  2  Quartzite,  Rutland,  Vermont. 

3  Granite,  Westerly,  Rhode  Island.  4  Limestone, 
Oak  Hall  Quarry,  Pennsylvania.  6  Fused  quartz. 


Vd,  Void  Rotio 


Figure  9.  G*  vs  void  ratio.  1  Granite,  Stone  Moun¬ 
tain,  Georgia.  2  Quartzite,  Rutland,  Vermont. 

3  Granite,  Westerly,  Rhode  Island.  4  Limestone, 
Oak  Hall  Quarry,  Pennsylvania.  6  Fused  quartz. 
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Q  values  for  frozen  soils  were  found  in  the  range  of  10-100.  They  are  bounded  below  by 
polycrystalline  ice  and  above  by  solid  rock,  such  as  quartzite  (Q  =  250)  and  granite  (Q  =  70) 
(Volarovich  and  Gurvitch,  1957). 

Although  many  theories  concerning  the  attenuation  of  stress  waves  in  earth  materials  have 
been  proposed,  none  of  them  is  definitive.  The  attenuation  mechanism  is  difficult  to  pin  down 
even  in  the  laboratory,  where  measurements  can  be  made  as  a  function  of  frequency,  temperature, 
pressure,  purity,  grain  size,  annealing  history,  etc.  In  the  present  work  we  used  viscosity  to 
describe  anelastic  properties  of  frozen  soils.  When  applied  to  solidn,  this  term  usually  means  that 
stress  relief  and  deformation  occur  by  some  poorly  understood  process,  which  may  be  a  combination 
of  several  types  of  processes,  and  may  result  in  linear  or  nonlinear  internal  friction  of  complicated 
(or  unknown)  frequency  dependence  (Jackson  and  Anderson,  1970).  Further  efforts  are  required  to 
obtain  complete  understanding  of  attenuation  in  frozen  soils. 


Determination  of  Linear  Viscoelastic  Constitutive  Equations  from  Measured  Complex  Modulus 

The  general  constitutive  equations  of  a  linear  viscoelastic  solid  are  given  as  (Eringen,  1967) 

R<rrAl  +  S<7kl  =  P<rAl  +  2®fkl  ^ 

where  <7^,  are  stress  and  strain  tensors  respectively,  and  P,  Q,  R  and  S  are  linear  operators 
defined  as. 


pa  J1  \  ** 

P  =  AQ  +  2  A,  - 
i=  1  dtl 


n2  # 

<?  =  Po  4  1  Pi  — 
«=1  <9t* 


R  =  ar 


+  Z.  a. 


«=1  dP 


*4 

S  =  Po  +  1  Pi  — 
i  =  l  dP 


(10a) 


(10b) 


(10c) 


(10d) 


where  Ajt  gt,  and  Px  are  constants. 

Torsional  mode 

Theory.  We  deal  only  with  circular  cylindrical  specimens  and  use  cylindrical  coordinates  for 
convenience.  Let  the  coordinates  be  r,  6  and  z,  with  z  being  the  axis  of  the  cylindrical  specimen, 
and  let  the  corresponding  displacements  be  uf,  uq  and  uz.  In  the  propagation  of  torsional  waves,  no 
longitudinal  or  lateral  displacement  is  to  be  expected  and  the  motion  is  symmetric  about  the  axis  of 
the  cylinder.  Therefore,  uf  and  uz  must  both  vanish  and  we  need  to  consider  only  the  wave  equation 
for  it q.  If  the  torsional  stress  applied  to  the  specimen  varies  sinusoidally  with  time  and  the  strain 
thus  induced  also  varies  sinusoidally  but  with  a  phase  difference,  then  we  may  write  the  wave  equa¬ 
tion  for  viscoelastic  materials  as: 
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d2ud  dzud 

=  <11) 

where  G  is  a  complex  shear  modulus,  p  is  the  density  of  the  specimen,  and  t  is  time. 

The  wave  equation  is  also  written  in  terms  of  the  present  linear  viscoelastic  model  as: 

dhe  da0z 


dt‘ 


dz 


SaQz  =  Q 


du 


0 


dz 


(12a) 


(12b) 


Among  many  alternatives  we  selected  the  following  four*parameter  model  for  torsional  vibration: 

(13a) 


Q  d  dz 

Q  ■  "i  s  +  * 


dtd 


dz_ 

dt2 


(13b) 


If  Uq  and  oqz  are  harmonic  with  an  angular  frequency  w,  then  the  complex  shear  modulus  G  is  given 
as 


C  =  Cj  +  i  G2 


where 


C1  - 


"  ^2^0  "  ^2w2)1 


(14) 


(15a) 


w[/ij(/30  -  /32&>2)  +  p2a>2] 


u3 

h  =  (fi0  -  ^2w2)2  + 


The  phase  velocity  V,n  and  the  group  velocity  V are  obtained  as  follows: 


%  = 


tg 


K 

\  1  (Gl\  1 

COS  —  -  —  (O  I -  I  COS  —  -  - 

2  2  \G*  /  2  2 


1  s  .  t 
-  ~  w  Sln  ~ 


(15b) 

(15c) 


(16a) 


(16b) 
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where 


G*  =  j  G 


(17a) 


G 


♦ 

OJ 


(17b) 


(17c) 


(17d) 


The  values  of  constants  /j1(  /j2,  /30  and  /Sg  are  determined  from  the  four  observed  variables,  namely 
the  complex  modulus  and  the  loss  factor  Lg(=  tan  <5t),  at  any  two  distinct  frequencies.  Since  it  is 
not  generally  expected  that  the  model  satisfies  these  four  conditions  exactly,  the  following  scheme 
is  used.  Suppose  the  observed  modulus  and  the  loss  factor  at  two  distinct  frequencies  are  G^\ 
G<2>,  L(sl),  and  L*2)  respectively;  then  find  /jj  and  nz  which  minimize  the  derivation  V v  defined  as; 

l(G<l>  -  G<1>)2  +  (C<2>  -  G<2>)2] 

v  = - - - 52 -  (18) 

{G(1)2  +  C(2)2) 


under  the  constraints 


sm  s 


(19a) 


'■iS  ‘  (19») 

where  Gm  and  Lgm  are  the  complex  modulus  and  the  loss  factor  predicted  by  the  model. 

Results  and  discussion.  The  four  parameters  n  j,  nz,  Pq  and  /3g  were  computed  as  follows.  We 
used  the  observed  modulus  and  the  loss  factors  at  1.00  kHz  and  2.00  kHz  to  compare  the  constitutive 
equations  of  different  frozen  soils  around  1.5  kHz  frequency.  Since  eq  19a  and  19b  are  linear  in 
terms  of  PQ  and  we  solved  these  equations  to  substitute  and  nz  for  PQ  and  Pz  in  eq  18.  Now 
V i  is  a  function  of  /jt  and  nz  only;  that  is, 

vt  =  vM>Pz)-  (2°) 

When  V{  is  a  minimum,  we  have 

°  (21a) 

dp  x 

dVl 

G(nx.  H2)  =  0. 

6tiz 


(21b) 
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Suppose  that  is  an  approximate  solution  and  let  8p^  and  8p.^  be  coiTections 

which  we  shall  determine.  Expanding  F  and  G  in  a  Taylor  series  and  truncating  after  the  first- 
order,  we  get: 


F(n\l\  /41})  +  8p[1) 


©<#»?>•  #41>) + 


4”)  »r<A 1>-41>) 


ty  1 

ac^1).  41}) 

i 


dFs 


=  0 


<9g(41>,  41)) 

541) - 1 — —  =  0. 


(22a) 


(22b) 


This  linear  system  in  8p^  and  8^  gives  the  next  approximation,  ^  and  as 


li[Z)  =  41}  +  s41}  (23a) 

42)  =  41}  +  S41}-  (23b) 


The  above  iterative  procedure  is  repeated  until  we  obtain  p  j  and  which  minimize  l',.  The 

results  of  the  computation  are  shown  in  Table  II. 

CGS  units  were  used  in  the  computation.  For  instance,  stress  was  expressed  in  dynes/cm3. 

The  computation  was  quite  satisfactory.  An  absolute  minimum  value  of  Vt  exists  for  all  specimens 
computed  and  convergence  is  rapid.  In  order  to  indicate  the  degree  of  approximation  or  error  we 
listed  V'f  in  the  table.  The  values  of  v'f  vary  in  the  range  of  10‘2  to  10*3,  which  is  considered 
satisfactory.  In  the  general  constitutive  relation,  eq  9,  if  all  other  constants  except  Aj,  p1  and  /3j 
are  vanishing  and  is  equal  to  unity  as  defined,  eq  9  reduces  to  the  elastic  constitutive  equation, 
or  Hooke’s  law  in  terms  of  an  incremental  strain  resulting  in  an  incremental  stress,  where  Aj  and  pl 
are  Lamd  constants.  Therefore,  p2,  and  /3g  indicate  a  degree  of  deviation  from  elastic  solid. 

For  ice-saturated  or  almost  saturated  non-plastic  frozen  soils  a  strong  correlation  was  found  between 
modulus  and  void  ratio  as  described  before.  We  plotted  four  parameters,  pv  p2,  /3Q  and  /92,  versus 
void  ratio  in  order  to  find  out  some  general  trends  (Fig.  10).  It  appears  that  a  correlation  exists 
between  either  p1  or  p2  and  void  ratio,  but  neither  /SQ  nor  /8g  are  affected  by  variations  of  void  ratio. 

Longitudinal  mode 

Theory,  If  the  wavelength  of  the  harmonic  wave  is  long  in  comparison  with  the  diameter  of  the 
specimen,  the  equation  of  motion  for  longitudinal  vibration  in  Cartesian  coordinates  is  written  as: 

d2u  dzu 

E - I  =  p - 1  (24) 

dxz  dtz 


where 


E  =  complex  Young’s  modulus 
p  =  density  of  the  medium 
t  =  time 

ux  =  displacement  along  x  coordinate. 
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In  this  case  axx  is  the  applied  stress  and  the  other  five  components  of  stress  are  zero.  The 


first  three  equations  in  eq  9  thus  become: 

(ft  +  S)axx  =  P  A  +  2  Qexx  (25a) 

R  CTXX  -  P  A  +  2Q<yy  (25b) 

Raxx  =  PA  +  2  Q(zz.  (25c) 

Eliminating  <yy  and  ezz  from  eq  25.  we  obtain  the  following  equation: 

[Q(R  +  S)  +  PS]  o xx  =  (3 PQ  +  2QZ)(XX.  (26) 


The  operators  Q  and  S  have  been  determined  and  the  operators  P  and  ft  are  now  to  be  determined 
from  the  longitudinal  vibration  tests.  The  determination  of  P  and  ft  is  more  elaborate  than  the 
determination  of  Q  and  S,  since  restrictions  have  already  been  imposed  on  the  former  in  relation  to 
the  latter. 

We  intend  to  introduce  four  new  parameters  to  define  P  and  Q.  There  are  several  ways  of 
choosing  such  parameters.  In  the  general  constitutive  relation,  eq  9,  if  all  other  constants  except 
Aj,  and  /3j  are  vanishing  and  /3j  is  equal  to  unity  as  already  defined,  eq  9  reduces  to  the  elastic 
constitutive  equation,  or  Hooke's  law  in  terms  of  an  incremental  strain  resulting  in  an  incremental 
stress,  where  Aj  and  (i j  are  Lam£  constants.  It  is  anticipated  that  the  dynamic  properties  of  frozen 
soils  do  not  deviate  markedly  from  those  of  an  elastic  solid.  Thus  the  parameter  Aj  is  expected  to 
play  an  important  role  in  the  operator  P.  We  choose  two  models  defined  by  two  four-parameter  groups, 
(Aq,  Aj,  Oj,  ag)  and  (Aj,  Ag,  aj,  ag).  Now  we  have 
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P 


= 


dZ 

di2 


R  - 


d2 
dt 2 


(27a) 


(27b) 


where  A2  =  0  for  Model  1  and  AQ  =  0  for  Model  2. 

When  ctxx  and  r  are  harmonic  with  an  angular  frequency  w,  then  the  complex  Young’s  modulus 
E  is  given  as. 


£  +  iAfg 

L  |  +  i  L  2 

where 

Lj  =  +  Il2a>3a2  ~  +  ^2^0  +  ^1^1  + 

+■  Ag^o)  +  to^gg/^  +  Ag/)g)  f  A0(/30  —  fipco2)/ co 
^2  =  ~  ^2a>‘“a  1  "  ^1£l»2«2  +  @0^1  f  ^1^0  +  ^0^1  - 

-  tL»2(^ij^2  +  ^2^1  +  + 

Mj  =  -tuiSAj/ij  +  2/ij  +  3Aq//2)  +  w3(  3Ag/x2  +  2/i|)  (29c) 

M2  =  3A0/ij  -  w^Ajfig  +  3A2/Xj  +  4njfig).  (29d) 

Knowing  the  values  of  /ij,  fig,  /^o  an<*  ^2’  we  determine  two  groups  of  four  unknown  constants, 

(Aq,  Aj,  Oj,  ag)  and  (AJt  A2,  ay  a2)  from  the  observed  complex  modulus  and  the  observed  loss  factor 
at  any  two  distinct  frequencies  by  a  method  similar  to  that  used  for  the  torsional  mode.  Finally  we 
select  one  of  the  groups,  which  minimizes  the  error  Fj,  defined  as: 

l( E(1)  -  £(1))2  +  (£(2)  -  E(2))21 

F,  - - JH -  -  (30) 

(E(i)2  ,  £(2)2} 


(28) 


(29a) 


(29b) 


where  E('\  E^’  are  defined  as  in  torsional  vibration.  The  phase  velocity  F1[}  and  the  group  velocity 
Fj  are  given  as  for  the  torsional  mode: 
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where 

E_ 2 
E1 


E{  =  Re  (£) 

£2  =  Im(E). 

Results  and  discussion.  The  actual  computation  of  four  parameters,  either  (A0,  Aj,  a1(  a2)  or 
(Aj,  Ag,  aj,  a2),  was  made  in  the  same  way  as  in  the  torsional  mode.  We  used  the  observed  modulus 
and  the  loss  factors  at  1.00  kHz  and  2.00  kHz.  It  turns  out  that  Model  1  always  gives  a  better 
approximation  than  Model  2  for  all  ice-saturated  frozen  soils  examined.  The  results  of  the  computa¬ 
tion  are  shown  in  Table  III.  The  degree  of  approximation  is  not  satisfactory.  The  values  of  V^ 
vary  in  the  range  of  10‘J  -  10'\  which  is  much  larger  than  v'f.  It  might  be  possible  to  obtain  a 
better  or  well  balanced  approximation  for  both  torsional  and  longitudinal  modes  by  selecting  eight 
parameters  in  a  different  manner.  Also  one  could  improve  accuracy  by  introducing  more  parameters. 

The  most  commonly  used  linear  viscoelastic  models  for  earth  media  have  been  the  Maxwell 
or  Voigt  models  with  only  two  parameters  and  very  little  attention  has  been  directed  toward 
examination  of  the  complete  constitutive  equation.  Although  the  dynamic  behavior  of  frozen  soils 
is  less  complex  than  that  of  unfrozen  soils,  the  degree  of  deviation  from  perfect  elasticity  is 
surprisingly  large.  It  is  felt  that  further  efforts  should  be  made  to  investigate  the  constitutive 
equation  of  frozen  earth  materials  according  to  the  theory  of  continua. 


E*  =  \E\ 

<5j  --  tan-1 


r*  SE* 

£“  "  "ST 

dS , 

o  1 
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Table  III.  Values  of  constants  for  longitudinal  mode. 


( '  SS  1053  Manchester  Silt 


Ao  0.1 

A.  -4.2B  X  10" 

a.  7.05  x  10  "’ 

*.  6.50  x  t0'( 

Vx  I.53  x  10'* 

v/4>*  1 .41  X  10’ 

1.0  _ 
-4.28  x  10* 

7.60  x  10'*’ 

6.52  XlO" 

1.53  x  10** 

1.41  x  10*' 

5.0  n 
-4.71  x  10 

7.82  x  10' 
7.89  x  10" 

1.82  x  10-* 
I.76  x  10*' 

Id Ottawa  Sand 

.At  0.1 

>,  1.50  x  10'* 

A.  4.80  X  10*“ 
v  -1.92  x  10"' 

I.87  x  10'* 

Vi'A  1.18  X  10*' 

1.0 

1.41  x  10'  * 

4.52  x  10" 

3.21  x  10  * 

2.39  x  10*‘ 

1.10  x  10*'' 

1.33  x  10  * 
4.21  x  10'* 
5.65  x  10' 

3.83  x  10“ 
1.02  x  10"' 

Si\  1027  Fairbanks  Silt 

v'e  0.1 

A.  -1.44  x  10'1 

A.  4.37  x  IO" 
at.  1.77  x  10 

4.40  x  10 ** 

6.64  x  IO2 

1.0  ^ 

-7.13  x  10 

3.73  x  10" 
-7.04  x  10'J 

1.28  x  10*' 

3.78  x  10'2 

5*0 

-1.01  x  10' 
9.13  X  10" 
5.58  x  10"' 

2.35  x  10  ' 

8.37  x  10  ' 

SN  1024  Polycrystalline 

Ice 

<*■*  0.1 

A.  4.26  x  10,+ 

A.  4.58  X  10'° 

X.  -6.75  x  10'* 

•X.  1.36  X  10'' 
v/  1.32  x  10-' 

i.o  /v. 

1.91  X  10'* 
4.56  x  10" 
-5.37  x  10*' 

3.70  x  10*‘ 

8.35  x  1C' 

5.0 

4.60  x  10'* 
4.79  x  10" 
-2.35  x  IO" 
9..41  x  10  J 
1.22  x  10  ' 

SN  1023  Polycrystalline 

Ice 

An  o.l 

A.  3.82  x  10"* 

A.  4.69  x  10'“ 

*>  -5.73  x  10  ' 
of.  3.31  X  10"‘ 

Vi*  6.65  xlO*> 

uo  w 

4.38  x  10' 

4.69  X  10'“ 
-7.01  x  io*; 

7.55  X  10* 

1.02  x  10*' 

3.05  xl0'*o 
4.34  x  IO10 

-1.39  x  io;t 

5.44  x  10  , 
4.67  x  10  "' 
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ID.  THE  USE  OF  FREE  OSCILLATIONS  TO  MEASURE  THE  ELASTIC 
PROPERTIES  OF  MATERIALS 

by 

M.  Smith,  R.  Martin,  and  Y.  Nakano 


Introduction 

This  section  discusses  our  efforts  to  develop,  and  apply,  experimental  techniques  for  the 
measurement  of  the  elastic  and  slightly  anelastic  properties  of  materials  based  on  the  free 
resonances  of  layered  elastic  spheres.  We  initially  undertook  this  effort  because  such  spheres  are 
the  only  bounded  bodies,  known  to  us,  for  which  we  could  develop  exact  analytic  solutions.  Con¬ 
sequently  we  are  not,  a  priori,  restricted  to  considering  only  high  or  low  frequency  approximations 
or  time-limited  response. 

Although  some  classic  studies  of  the  dynamic  properties  of  elastic  spheres  have  long  been 
available  (see  Love,  1944),  the  results  were  generally  devoid  of  practical  significance  until  the 
Earth’s  free  oscillations  were  observed  after  the  Chilean  earthquake  of  1960  (Bullen,  1963).  The 
consequent  attention  by  seismologists  resulted  in  a  well-developed  literature  from  which  many  of 
our  references  are  drawn.  We  differ  slightly  in  that  we  restrict  our  attention  to  spheres  com¬ 
posed  of  discrete  layers  and  we  also  neglect  the  effects  of  self -gravitation,  rotation,  an  initial 
stress  state,  and  ellipticity  (see  Dahlen,  1968). 

The  following  three  sect;ons  and  the  appendices  are  theoretical,  with  some  numerical  examples. 
So  far  as  we  are  aware,  the  particular  development  given  here,  that  is,  a  layered  non-gravitating 
sphere,  has  not  been  published,  of  a  piece,  elsewhere.  It  is,  however,  a  “standard,”  albeit 
complicated,  problem.  We  believed  that  its  detailed  solution  had  to  be  explicitly  laid  down  before 
progressing  into  experiment. 

The  last  sections  deal  with  the  results  of  pilot  experimentation.  We  believe  the  results  indi¬ 
cate  that  the  practical  difficulties  associated  with  this  method  are  being  mastered  and  that  the 
technique  is  a  viable  one. 


Elastic  Displacement  Solutions  in  Spherical  Coordinates 

We  consider  a  volume  of  space  filled  with  an  isotropic,  homogeneous,  linearly  elastic  medium 
having  Lamd  constants  A  and  p,  and  density  p.  We  assume  the  medium  to  be  free  of  gravitation  and 
other  body  forces,  but  allow  the  existence  of  one  of  more  surfaces  across  which  tractions  may  be 
applied. 

Let  ube  the  displacement  field  specifying  the  motion  of  each  particle  from  its  unique  rest 
position.  We  assume  c  to  be  a  first  order  infinitesimal  and  do  not,  therefore,  have  to  distinguish 
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between  Eulerian  and  Lagrangian  coordinate  systems.  Let  T  be  the  elastic  stress  tensor.  If  we 
assume  that  u  =  0  corresponds  with  the  unstressed  state  of  the  medium,  T  is  given  by  (Fung,  1965) 

T  =  A(  V  •  u )/  +  /i(Vu  +  uV),  (1) 

where  l  is  the  identity  tensor.  Vu  is  the  gradient  tensor  of  u,  and  uVis  its  transpose. 

The  conservation  of  linear  momentum  leads  immediately  to  the  equation  of  motion, 


pdf  u  =  V  •  T.  (2) 

Equation  1  and  some  standard  vector  calculus  identities  convert  eq  2  to 

pdf  u  =  (A  +  2/i)V(V  •  u)  -  /jV  x  V  x  u.  (3) 

We  choose  to  represent  u  by 

u  =  7 U  +  -7  x  VjW  (4) 

where  U,  V,  and  W  are  scalar  fields  and  Vj  is  defined  by 

Vi  =  03q  +  <p  (sin  6)~ 1 3<j> .  (5) 


7is  a  unit  vector  directed  away  from  the  origin,  0  and  cf>  are  the  colalitude  and  longitude,  and  0  and 
are  their  respective  unit  vectors.  Vj  is  the  gradient  operator  on  the  surface  of  a  sphere  of  unit 
radius.  It  is  related  to  the  three-dimensional  gradient  by 

V  =  7<Jr  -  r_1  Vj.  (6) 

After  some  algebra,  we  can  show  (Backus,  1967)  that 
V(V-  u)  =7dr  j(<9f  +  E)t/  +  + 

+  +  r-8vfFj>,  (7) 


and 


V  x  V  x  u  =  r 


tr  2(9r(rVjF)  -  r  +  Vj  lr 


-2„2 , 


.-1 


dtU  - 


■1j2 


df(rV)  I 


+  7  x  Vj  If  2V?W  +  r~ldf  Wi. 


(8) 


We  insert  eq  4,  7  and  8  into  eq  3.  We  now  appeal  to  the  uniqueness  of  the  representation  4  (Backus, 
1967)  to  yield  the  three  coupled  partial  differential  equations 
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pS^U  .  (A  *  2|0dr/(at  4  r-)  V  ,  r-‘vf vj-  - 

-  ^  U,(rvfl')  - 

r * 

p0?F  =  (A  +  2n)I(rldt  + 


-  £  [dTu  -  a2(rF)l,  (9b) 

r  r  r 

and 

pafw  -  fiafwo  +  rVfm.  (9c) 

We  note  that  both  V  and  W  may  be  augmented  by  ary  constant  without  affecting  a  (see  eq  4).  There¬ 
fore,  we  may  expect  these  two  scalar  fields  to  be  determined  only  to  within  an  additive  constant. 

We  also  note  that  eq  9a  and  9b  both  involve  U  and  V  but  not  W,  while  eq  9c  involves  W  but  neither 
of  V  and  V. 

The  manner  in  which  one  elects  to  solve  eq  9a,  9b  and  9c  (plus  whatever  boundary  conditions 
appertain)  depends  upon  the  intended  application  of  the  solution.  For  our  purposes,  we  wish  to  find 
a  set  of  complete,  linearly  independent  vector  fields  if1,  f2,  . . .  I  each  of  which  satisfies  eq  3.  If  the 
set  is  complete,  all  possible  solutions  to  eq  3  may  then  be  expressed  as  some  linear  combination  of 
the  members  of  the  set,  the  coefficients  used  in  the  expansion  being  determined  by  boundary  and 
initial  conditions. 

In  pursuit  of  this,  we  Fourier  transform  eq  9a-9c,  going  from  time  t  to  angular  frequency  w.  We 
do  not  introduce  a  distinct  symbol  for  Fourier  transforms  since  it  will  be  clear  from  the  context 
whether  a  symbol  refers  to  the  transformed  or  untransformed  variable.  The  result,  of  Fourier  trans¬ 
formation  is  to  replace  d2  by  -co~. 

To  transform  the  resultant  trio  of  equations  from  partial  to  ordinary  differential  equations  we 
introduce  a  surface  spherical  harmonic  expansion  of  U,  V  and  W.  For  l  >  0  and  -I  <  m  <  I,  we 
define  (Hill,  1953) 


vfl/ 1  ,  (9a) 

1)U  .  r-V*'}  - 


Y?  (0,  <f>)  =  (-l)m  |-J-t  1  jmj|!  1 Pf  (cos  0)e,m^  (10) 

1  I  4ff  (1  +  lm|)!J  1 

where  P™  is  the  Associated  Legendre  Function  given  by 

P™(x)  =  l1  ~  x2 
2>ll 

If  Sj  is  the  surface  of  a  sphere  o’  unit  radius  centered  on  the  origin,  the  Y “  are  orthonormal  in  the 
sense 

/  Yf  (0,  &  Y%  (0,  J)  sin  0d0d<£  =  SlAS  (12) 

sl 

whet e  the  bar  indicates  complex  conjugation  i  id  5,j  is  the  Kronecler  delta. 


)(|w|/2) 


a/Hrnl 


(x2  - 


l)1. 


(11) 
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The  Y (0,  (fj)  form  a  complete  set  and,  if  we  assume  U,  V  and  W  are  sufficiently  regular,  we 
may  expand  them  as 

U  (r,  6,  <f>,  co)  =  I  1  U?  (r,  u)  Y?  (0,  <f>),  (13a) 

1=0  m=-l 

V(r.  6,  <t>,  co)  =  1  2  V ®  (r,  w)  T®  (0,  <f>),  (13b) 

1-1  m=-l 

and 

W(i,  0,  4>,  oj)  =  1  I  W®  (r,  u>)  Y?  ( 0 ,  <f>).  (13c) 

1=1  m=— 1 


The  1  =  0  terms  have  been  omitted  from  eq  13b  and  13c  since  inspection  of  eq  4  reveals  that  these 
terms  do  not  contribute  to  the  displacement  field. 

We  now  insert  eq  13a*13c  into  the  transformed  equations  9a-9c,  and  make  use  of  the  relation 
=  -  HI  +  1)T® .  ,  (14) 


The  resultant  expressions  are  then  multiplied  by  a  particular  K®  and  integrated  over  the  surface  of 
a  sphere  of  radius  r.  Application  of  eq  12  leads  immediately  to 


poj2l/™ 


=  (A  +  2fi)dt  |(<9r  +  £){/™  -  iLlJl 


and 


+  JL  |i(I  +  l)at(rK®)  -  l{l  +  Dt/fl 

r2 

- 1*  *  **> {(*>,  *  •'"}  - 


ii  I 0U®  -  dfirt'®)!, 

r  I  1  I  1 


-  =  JL  |^(rW®)  -  1(1  L  1}  W®  J  . 


(15a) 


(15b) 


(15c) 


We  note  that  none  of  eq  15a-15c  explicitly  involves  m. 

The  set  15a- 15c  can  be  solved  by  any  of  several  standard  techniques.  The  method  used  here 
is  detailed  in  Appendix  A.  The  results  are 
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'tyjOr r)  —  +  —  j,(yr) 

» 

- r~1<9r[r/1(yr)] 


5ry,  Clcr) 
yx  (kr) 


HI  +  l) 


My') 


■^[r /,  (yr)] 


for  1  >  0; 


Uo  =  Wkr)1 


ft} 


and 


Vn  =  0; 


ft} 


tt'1  =  lr/j(yr)  ry^yr)!  <  V  for  /  >  0; 


Wn  =  0; 


where 


k  =  — 


Cl) 


(A  +  2p)  Vp  * 
P 


anu 


CO  CO 


(16) 

(17a) 

(17b) 

(18a) 

(18b) 

(19) 

(20) 


Vp  and  Vg  are,  respectively,  the  compressional  and  shear  velocities  of  the  medium.  Each  of  Ay, 

By . Fy  represents  a  spherical  harmonic  of  degree  1.  That  is,  each  one  is  some  linear  combina¬ 
tion  of  the  21  +  1  functions  . y®,  . . y{l  but  we  cannot  specify,  without  considering  a 

particular  problem,  what  linear  combination  each  is.  To  be  specific,  we  may  express  them  as 

Ay(0,<f>)  =  i  A®y®(0,<£),  (20a) 

m=-l 


and  similarly  for  By,  Cy,  Dy,  £t  and  Fy.  The  A, . . .  are  not  constants  but  the  4®  . . .  are.  Equations 
16-18  thus  contain  6 (21  +  1)  presently  undetermined  constants.  (When  1  =  0,  there  are  only  two 
constants.) 

We  believe  that  the  manner  ir.  which  we  have  dealt  with  the  spherical  harmonics  merits  further 
elaboration.  We  will  use  eq  10  as  an  example.  We  could  retain  the  superscript  m  and  have 


W®  =  \B?rjy(yr)  +  F®rMyr)iy®(0,<£). 


(21) 
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Let  us  then  define 

wl  =  {  W®. 

m=— 1 

Combining  these,  we  have 

W,  =  r;,(yr)  I  Ef  Yf  (6,  *)  +  ry,  (yr)  1  F®y®(0,  <£) 

m=-l  m=-l 

where  we  were  able  to  regroup  the  sums  because  each  of  the  two  independent  solutions  to  eq  15c, 
namely  rj j  (yr)  and  ryj  (yr),  is  independent  of  m.  This  last  expression  is  simply  eq  18a. 

Any  linear  combination  of  spherical  harmonics  of  the  form 

Hx  =  i  a®K®(0,<£)  (22) 

m=-l 

is  itself  a  spherical  harmonic  and  satisfies 

vff/j  -  1(1  +  l)Hj.  (23) 

The  set  IF^,  . . .,  F° . y[\  serves  only  as  an  orthonormal  basis  for  the  21  +  1  dimensional 

space  of  spherical  harmonics  satisfying  eq  22,  and  has  no  inherent  significance.  At  this  point  in 
our  development  we  can  only  know,  then,  that  each  set  of  radial  functions  in  eq  16-18  will  be 
multiplied  by  some  spherical  harmonic,  H j.  We  cannot  know  H j’s  precise  form;  the  a®  in  eq  22  will 
be  determined  by  boundary  and  initial  conditions. 


Fres  Oscillations  of  a  Layered  Elastic  Sphere 

We  consider  a  sphere  divided  into  N  concentric  spherical  shells.  We  number  the  shells  from 
the  center  outwards  and  let  be  the  outermost  radius  of  the  ith  shell.  Then  rN  is  the  radius  of  the 
sphere.  Let  r0  equal  zero.  We  suppose  that  the  ith  shell  is  composed  of  an  elastic,  homogeneous, 
isotropic  medium  of  density  px  and  having  Lamd  parameters  and  px.  These  parameters  define  the 
compressional  velocity  Fpj  and  the  shear  velocity  Fsi. 

We  assume  that  the  surface  r  =  rN  is  free  from  all  tractions,  and  that  no  body  forces  (such  as 
gravitation)  are  present.  We  wish  to  know  for  what  angular  frequencies  a>  there  exists  a  non-trivial 
displacement  u(r)  elw£,  such  that  the  surface  is  free  of  traction  and  all  internal  boundary  conditions 
(discussed  below)  are  fulfilled.  We  shall  label  such  angular  frequencies  eigenlrequencies  and  their 
associated  displacements  eigenfunctions.  We  refer  to  such  traction-free  motions  as  free  oscillations. 

Let  u^V)  eiwt  be  the  displacement  field  in  the  layer.  From  the  displacement,  we  can  compute 
a  stress  tensor,  T(i\  by  eq  1.  Equation  2,  namely, 

pdf  0  =  V-  T  (2) 

must  hold  everywhere  in  the  medium,  since  it  expresses  only  the  conservation  of  momentum  and  is 
not  dependent  upon  such  assumptions  as  isotropy,  homogeneity,  etc.  In  particular,  eq  2  is  valid 
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in  a  small  “Gaussian  pillbox”  which  encompasses  a  portion  of  the  surface  r  =  rjt  the  boundary 
between  the  itb  and  (i  +  l)tft  shells.  Let  v  denote  the  volume  enclosed  between  the  radii  ri  -  Sr 
and  rt  +  Sr  and  by  some  range  of  the  coordinates  6  and  <f>.  Then 


f  pd*  udv  =  /  y.  Tdv.  (24) 

v  V 

If  1  is  the  surface  of  v.  Gauss's  theorem  leads  to 

/ pdfudu  =  J  T  •  lido 
v  1 

where  n  is  the  unit  outward  normal  on  2 .  If  we  let  Sr  go  to  zero,  then  the  volume  of  v  also  vanishes 
and,  if  p  u  remains  bounded,  the  left-hand  integral  goes  to  zero.  The  right-hand  side  becomes 


f  |T (j+1>  -  T(i_1)  |  .Tda  =  0  atr  =  r,.  (25) 

1 

Since  this  result  is  independent  of  the  details  of  the  shape  of  1  we  conclude  that  the  quantity  T  •  T 
must  be  everywhere  continuous,  and  in  particular  across  boundaries. 

To  condition  25  we  add  one  expressing  our  intuition  of  the  behavior  of  elastic  materials.  If 
both  the  ith  and  (i  +  \)th  shells  are  solid,  we  require  that  the  displacement  n  be  continuous.  We 
refer  to  interfaces  at  which  this  is  true  as  being  “welded.”  If,  however,  one  or  both  shells  are 
fluid  (that  is,  p  =  0),  we  require  only  the  radial  component  of  displacement  to  be  continuous.  In  the 
latter  case,  we  allow  the  boundary  to  slip  but  in  neither  case  do  we  allow  holes  to  open  or  matter  to 
interpenetrate  itself. 

The  quantity  T  •  Tis  a  vector  and  represents  the  traction  (force)  acting  on  a  surface  normal  to 
77  In  a  fashion  identical  to  eq  4,  we  may  represent  it  as 

T  .  T  =  TP  +  VtQ  -  T  x  VjR,  (26) 

where  P,  Q  and  R  are  scalar  fields.  Equation  25  states  that  P,  Q  and  R  are  continuous  across  an 
interface.  The  stress-displacement  relations,  eq  1,  enable  us  to  relate  P,  Q  and  R  to  U,  V  and  W, 
the  scalar  representatives  for  u.  These  relations,  which  are  derived  in  Appendix  B,  are 

P  =  (A  +  2  n)drU  +  —  U  +  (27a) 

1  r  r 

Q  ■  l‘{~,  * ,3'  (•;)}•  (87b) 

and 

R  =  ^(7)-  <27c) 

The  import  of  the  boundary  conditions,  then,  is  that  P,  Q,  R  and  U  are  continuous  everywhere 
and  V  and  W  are  continuous  in  solid  domains. 
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We  now  have  six  scalar  fields  to  contend  with.  However,  an  examination  of  eq  16,  17,  18  and 
27  indicates  that  we  can  group  them  into  two  sets,  one  consisting  of  U,  V,  P  and  Q,  and  one 
consisting  of  5F  and  R.  These  two  sets  are  completely  independent;  they  do  not  interact  in  any 
way.  And  we  may,  without  loss  of  generality,  treat  them  separately. 

The  set  (U,  V,  P,  Q)  is  the  set  of  spheroidal  variables.  The  displacements  described  by  this 
set  give  rise  to  both  coi-.pressional  and  shear  strain,  and  are  interference  products  of  compressional 
and  vertically  polarized  shear  waves.  Among  other  things,  this  set  gives  rise  to  the  spherical  type 
of  Rayleigh  surface  waves. 

The  set  (W,  R)  is  the  set  of  toroidal  variables.  The  displacements  described  by  this  set  are 
orthogonal  to  t  and  produce  only  shear  strains.  They  are  the  interference  products  of  horizontally 
polarized  shear  waves.  Love  surface  waves  are  associated  with  this  set. 

From  this  point  forward  we  will  discuss  only  spheroidal  types  of  motion.  A  similar  development 
for  toroidal  oscillations  can  be  easily  formulated  since  the  toroidal  problem  is,  in  fact,  substantially 
simpler. 

For  any  specified  angular  frequency  of  oscillation  w,  the  set  ( U ,  V,  P,  Q)  can  be  expanded  in 
terms  of  spherical  harmonics  as 


U  (r,  6,  (fi,  C<j)  =  1 

1 

1 

U™  (r,  o)  Yf  (6,  4>) 

(13a) 

1=0 

m=-l 

V(r,  0,  <f>,  <o)  =  1 

1 

1 

Vf  (r.  <o)  Y?(6,  $) 

(13b) 

1=1 

m=-l 

P(r,  6,  4>,  o)  =  2 

1 

P™  (r,  w)  Yf{6,  0) 

(28a) 

1=0 

m=-l 

oc 

l 

Qf(r.  a.)  Yf(0.  Q. 

Q(r,  0,  <t>,  «)  =  1 

1 

(28b) 

1=1 

m=-l 

The  U J"  and  Fj"  are  given,  in  terms  of  a  set  of  coefficients,  by  the  analytic  solutions  16  and  17. 
P™  and  Q™  are  then  given  by  eq  27.  Our  problem  is  to  determine  those  frequencies,  o>,  for  which 
we  can  construct  U,  V,  P  and  Q  by  this  method,  for  each  layer,  such  that  all  boundary  and  interface 
conditions  are  met. 

However,  because  the  F™  are  an  orthogonal  set,  we  may  consider  the  above  problem  separately 
for  each  F™.  That  is,  given 


U(r,  0,  <f>,  to)  =  l/f(r,  w)  Y?(9,  <A),  (29) 

etc.,  for  what  angular  frequencies  «  can  we  satisfy  all  internal  and  external  boundary  conditions? 
Since  m,  as  discussed  earlier,  is  a  degenerate  index,  we  can  simplify  eq  29  to 

U  (r,  0,  <f>,  co)  =  (/j  (r,  o>)  Hl  (0,  tf),  (30a) 


V (r,  0,  4,u>)  =  Vx(j,u>)Hx{d,<!>), 


(30b) 
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P(r,  e,  0,  61)  =  P,(r,  to)  Hx  ( e ,  0),  (30c) 

and 

Q(r,  0,  0,  6>)  =  Qj  (r,  <*>)//,  (0.  0),  (30d) 

where  H }  is  given  by  eq  22  and  is  some  spherical  harmonic  of  degree  l.  The  details  of  H x  are  of 
no  particular  interest  since  the  eigenfrequencies  and  the  form  of  Ux,  . . .,  <?j  remain  fixed  no  matter 
what  H j  we  choose. 

Our  problem  now  is  to  find  the  eigenfrequencies  associated  with  spherical  harmonics  of  degree 

1.  If  we  wish  to  know  all  eigenfrequencies,  we  must  repeat  this  procedure  for  each  of  /  =  0,  1, 

2,  ... 

We  will  devise  a  constructive  algorithm  which  will  enable  us,  for  a  given  angular  frequency  o, 
to  construct  a  solution  from  the  center  outward  which  meets  all  boundary  conditions  save  one.  If 
the  last  condition  is  met,  is  then  an  eigenfrequency. 

We  consider  eq  16,  which  expressed  U j  and  Vx  as  a  linear  combination  of  four  independent 
functions  of  radius.  The  coefficients  Ax,  . . .,  Dj  were  taken  to  incorporate  all  spherical  harmonic 
content.  Equation  27  allows  us  to  extend  eq  16  to 


where  i  =  1,  ....  N  designates  the  layer  to  which  the  solution  is  appropriate.  H{  is  a  4  x  4  matrix 
constructed  from 

h3j  =  (A  +  2^)dthXi  +  -  A l{L±JlhSj,  (32a) 

and 


h4j  =  ftl  +  rdt(r~xh2j)  \  (32b) 

and  the  first  two  rows  of  Wj  are  taken  from  eq  16.  The  set  U|,  ...  d|!  serve  as  constants,  and  not 
spherical  harmonics.  For  the  remainder  of  this  development,  we  shall  omit  the  6  and  0  terms  for 
convenience. 

We  rewrite  eq  31  as 

S'  (r)  =  Hj(t)  .  C*1  (33) 

where  we  have  omitted  the  subscript  l.  The  vector  S'(r)  includes  both  the  stress  and  displacement 
terms. 

We  will  now  proceed  to  construct  a  solution  satisfying  all  internal  boundary  conditions.  In 
region  1  which  includes  the  point  r  =  0  we  can  a  priori  eliminate  those  solutions  which  go  as 
yj(kr)  and  yj(yr).  Therefore,  C1  has  the  form 
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where  e  t  and  are  four-dimensional  unit  coordinate  vectors  directed  along  the  first  and  second 
coordinate  axes.  In  the  first  region,  then,  we  have 

S^r)  H^r)  •  M1?!  +  B1^ '}.  (35) 

In  the  second  region,  we  have 

S2(r)  =  H2(r)  •  C2.  (36) 

Assuming  both  regions  to  be  solid,  the  boundary  conditions  require  that 

S2  (r  j)  =  S^r,).  (37) 

or 

_H2(rj)  •  C2  =  H1^)  .  M1^  +  BlT2\.  (38) 

Because  the  solution  composing  the  columns  of  B  is  linearly  independent,  matrix  theory  guarantees 
that  H  is  nonsingular.  Therefore,  we  may  express  C2  as 

C2  =  l//2(r1)]-1  .  tf^rj)  •  +  B1^!.  (39) 

An  alternative  form  for  eq  39  is 

C2  =  A1^2  +  B1^2  (40) 

where 

f2  =  L//2(r1)]— 1  •  tfVj)  •  e*i  (41) 

<2  =  lH2(r1)]~1  •  H1^)  -  T&.  (42) 

Equations  40,  41  and  42  suffice  to  specify  S2(r).  By  a  similar  procedure,  we  can  extend  the 
solution  from  the  ilh  to  the  (i  +  l)tA  shell.  The  appropriate  relations  are 


Si  +  1(r)  =  Hi+1(r)  • 

C'+1 

(43) 

Ci+1  Al(k+l  4 

BkCi+1 

(44) 

£i  +  1  =  [«i+1(r)r 1 

•  Bi  (r)  - 

(45) 

<Ti+1  =  [H,  +  1  (r)]  • 

H'  (r)  -  C- 

(46) 
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In  fact,  either  of  or  £i+1  can  be  computed  alone  by  starting  with  eq  41  or  42  and  simply 
applying  eq  45  or  46  as  many  times  as  is  necessary. 

We  suppoce  now  that  we  have  computed  (N  and  £N  where  N  is  the  number  of  shells.  The 
solution  in  this  shell  is  given  by 

8 N(r)  -  JiN (r)  •  \Al{"  +  (47) 

We  note  that  both  the  solution  associated  with  £N  [i.e.,  Si  (r)  when  B1  =  0]  and  the  solution 
associated  with  £N  separately  satisfy  all  of  the  internal  boundary  conditions.  If  w  is  an  eigen- 
frequency,  we  will  be  able  to  find  some  A1  and  B1,  not  both  zero,  such  that  the  last  two 
components  of  SN  (rN)  are  zero. 

A  straightforward  way  to  do  this  is  to  evaluate 

(a) SN  -  ifw(rN)  •  c;N,  (48a) 

and 

(b) Sw  =  Hw(rN)  •  (48b) 

If  Sg  (rN),  i.e.  P,  must  vanish,  A1  and  B1  must  satisfy 

A1(aSg )  =  -  B1  (bSg).  (49) 

We  may,  without  loss  of  generality,  impose  a  normalization  such  as 

(A1)2  +  (B1)2  =  1,  (50) 

which,  with  eq  49,  allows  us  to  compute  A1  and  B1.  Using  A1  and  B1,  we  then  compute  (rN),  i.e. 
Q,  and  examine  it  to  see  if  it  vanishes.  If  it  does,  <o  is  an  eigenfrequency;  if  it  does  not,  w  is  not 
an  eigenfrequency. 

If  oj  is  an  eigenfrequency  then  we  may  use  the  values  of  A1  and  B1  in  eq  43  and  44  to  compute 
the  eigenfunction  S i(r)  at  any  radius  in  the  system.  We  recall  that  S'  (r)  must  be  multiplied  by  some 
spherical  harmonic  of  degree  I  and  the  factor  eio)t  to  obtain  the  full  solution, 

Sj(r,  0,  </>,  t)  =  S{(r)Hj(0,  <f>)  eibJt 

where  we  have  reinstated  the  subscript  I.  We  emphasize  again  that  the  precise  nature  of  H j  depends 
upon  initial  conditions  and  is  not  relevant  here. 

This  method  requires  modification  when  either  /  =  0  or  I  £  0  but  one  or  more  shells  have  a 
vanishing  shear  modulus.  We  will  briefly  outline  the  form  these  modifications  take. 

When  1  =  0,  B1  vanishes  identically  and  only  the  solution  is  propagated.  The  matrix  H{  is 
collapsed  to  a  2  x  2  matrix  by  eliminating  those  solutions  with  arguments,  yr.  V  and  Q  both  vanish 
and  the  only  condition  at  r  =  rN  is  that  P  vanish.  A1  becomes  merely  a  scale  factor  and  may  be 
taken  as  equal  to  unity. 

In  a  fluid  region,  Q  vanishes  identically.  Across  a  solid/fluid  or  fluid/fluid  interface  V  may  be 
discontinuous  and  Q  is  continuous  and  zero.  If  we  are  propagating  upward  through  a  solid  and 
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Figure  11.  Loci,  in  the  (l,  l )  plane,  ol  the  (requencies  of  free  oscilla¬ 
tion  of  a  sphere  of  radius  10.16  cm  composed  of  material  having  a  density 
of  2.202  g/cm3,  a  compressional  velocity  of  5.968  x  10s  cm/sec,  and  a 
shear  velocity  of  3.764  x  10s  cm/sec.  Solid  lines  connect  modes  of  the 
same  overtone  number.  Dotted  line  indicates  the  asymptotic  slope  (dio/ 

<91  >  associated  with  Rayleigh  wave  propagation  in  an  infinite  half-space 
having  the  same  material  properties  as  the  sphere. 


encounter  a  fluid,  we  must  combine  the  £  and  £  solutions  to  yield  Q  =  0  at  the  interface.  This  re¬ 
quirement  determines  41  and  B1  and,  therefore,  C'  for  all  shells  up  to,  and  including  the  fluid 
region. 

Consequently,  in  a  fluid  region  we  have  only  one  solution.  In  crossing  from  a  fluid  to  a  solid 
region,  we  must  “start”  a  new  solution  having  some  non-zero  V,  but  for  which  U,  P  and  Q  vanish. 
We  may  always  find  some  £  which  yields  this  result. 

For  a  given  1,  we  shall  arrange  the  frequencies  of  free  oscillation  in  ascending  order,  as  0wlf 
jWj,  0wj,  ...  We  shall  designate  the  displacements,  as  a  function  of  as  0Uj,  jUj,  . . .  and  the 
four-vector  of  displacement  and  stress  by  0Sj,  jSj,  . . .,  etc.  The  "lowest”  mode,  for  a  given  1,  is 
referred  to  as  the  "fundamental”  mode  and  the  remainder  as  "overtones.” 


Figure  11  shows  the  loci,  in  the  (w,  l)  plane,  of  all  spheroidal  eigenfrequencies  lying  below 
100  kHz  for  a  homogeneous  sphere  with  10.16  cm  radius.  The  sphere  has  a  density  of  2.202  g/cm’, 
a  compressional  velocity  of  5.968  x  10s  cm/sec,  and  a  shear  velocity  of  3.764  x  10s  cm/sec. 

The  results  show  several  features  common  to  such  'alculations. 

The  lowest  fundamental  mode  is  0w2.  This,  almost  always,  is  the  case.  Secondly,  the  modes 
for  1  ^  2  tend  to  arrange  themselves  in  smoothly  varying  suites,  each  of  a  given  overtone,  or  radial, 
number.  These  are  kncwn  to  correspond  to  the  fundamental  and  higher-mode  Rayleigh  surface  waves. 
To  establish  the  connection,  we  observe  that  Y j(0,  <£)  has  the  value 


Y[(6,  <f>)  -  (-  l)i  I21  t ..ftl  eil^ 
1  y  4ff(2J)!  2l 


(51) 


as  can  be  seen  from  eq  10  and  11.  Y J  describes  motion  which  is  closely  confined  to  the  equator 
(0  -  ff/2)  and  which  behaves  as  waves  traveling  circumferentially  about  the  equatorial  zone.  The 
change  of  phase,  per  unit  of  distance  traveled  in  the  <£  direction,  is  i/a  and  is  therefore  the  mode’s 
surface  wave  number.  For  large  values  of  w,  we  may,  roughly,  expect  that  the  quantity 
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C/l  = 

dl 


(52) 


represents  a  group  velocity  and 


C*  wa 
1 


(53) 


represents  a  phase  velocity.  In  Figure  11  we  have  placed  a  dashed  line  whose  slope  is  given  by 
eq  52  when  U*  is  equal  to  the  group  velocity  of  Rayleigh  waves  propagating  along  a  half-space,  the 
properties  of  which  are  the  same  as  the  sphere’s.  We  regard  the  agreement  as  good. 

The  numerical  techniques  used  in  this,  and  subsequent,  calculations  are  discussed  in 
Appendix  C. 


The  Inverse  Problem 

In  the  preceding  two  sections  we  have  developed  the  “forward”  problem  for  a  layered  elastic 
sphere.  The  forward  problem  consists  of  the  generation  of  the  eigenvalue  spectrum  associated 
with  a  given  model.  In  this  section  we  consider  the  inverse  problem  of  utilizing  a  measured 
eigenvalue  spectrum  (which  will,  inevitably,  be  incomplete)  to  infer  the  properties  of  the  materials 
of  which  the  sphere  is  composed. 

Let  M  be  the  space  of  all  layered  elastic  sphsres  having  N  layers  delimited  by  the  points 

lr0,  rj . rN  I  where,  as  before,  rQ  =  0  and  rN  is  the  sphere’s  radius.  Then  all  models  in  M  have 

a  common  geometry  but  differ  in  the  elastic  properties  (including  density)  of  their  component  shells. 
M,  then,  is  i  space  of  dimension  3 N  (since  each  shell  has  three  distinct  properties)  and  we  may 
represent  a  g  ven  model  by  m,  a  vector  of  dimension  3 N.  We  further  limit  M  by  requiring  that  it 
encompass  on.  y  physically  realizable  models.  A  model  is  said  to  be  physically  realizable  if  each 
shell’s  properties  satisfy 

p  >  0,  (54  a) 

//  >  0,  (54b) 

and 

A>-f/i,  (54c) 

where  both  equalities  54b  aud  54c  are  not  simultaneously  true.  The  latter  two  constraints  merely 
express  the  condition  that  an  elastic  material  be  thermodynamically  stable  (Fung,  1965). 

Let  nWj  (m)  be  the  nth  overtone  of  the  spheroidal  mode  of  degree  l  associated  with  the  model  m. 
(Both  n  and  i  range  over  the  non-negative  integers.)  We  cannot  express  nw,  (m)  in  closed  analytic 
form  but  we  can,  through  techniques  previously  discussed,  generate  it  numerically.  We  can  now 
formulate  the  inverse  problem  in  the  followiug  manner. 

Let  (jjj,  i  =  1,  ....  if  be  observed  resonant  frequencies  associated  with  particular  modes  of 
oscillation.  We  wish  to  determine  a  model  m,  satisfying 


njWlj(m)  -  njO>[j 


i  =  1 . K . 


(55) 


y 
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As  Backus  and  Gilbert  (1967)  have  pointed  out,  we  do  not  know,  a  priori,  if  the  set  of  solutions 
to  eq  55  is  empty,  has  a  single  member,  or  is  a  subspace  of  M  of  one  or  more  dimensions. 

We  do  not  know  a  direct  procedure  for  solving  eq  55  for  one  or  more  models  m.  We  resort  here 
to  iterative  methods  for,  hopefully,  generating  successively  improved  approximations  for  in.  For 
convenience  we  rewrite  eq  55  as 

D{  (m)  D°  i  =  1 . K  (56) 

where  the  D°  are  data  and  the  D(  (m)  are  data  functions.  Dt  (m)  is  a  scalar-valued  function  whose 
domain  is  the  3Af-dimensional  vector  space  M  and  whose  value  is  the  angular  frequency  of  free 
oscillation  of  the  n[A  overtone  of  degree  /4.  Let  m0  be  some  model  which  we  beheve  to  lie  near  m. 
We  wish  to  find  some  perturbation,  5m,  in  m°,  such  that  m°  +  Sit  more  nearly  satisfies  eq  56. 

(m°  and  m°  +  5m  must  both  lie  in  M  but  5m  alone  need  not.)  We  wish  to  have 


D^m0  +  5m)  =  D?  i  =  1,  ...  if.  (57) 

Expanding  D^m)  in  a  Taylor  series  gives 

„  3 N  IdD , 

DAn°  +  Sm)  =  D.  (m°)  +  I  — 

1  1  f-'lAn. 

j=i  J 

Then,  to  first  order  in  |5m|,  we  wish  5m  to  satisfy 

dm. 

If  |5mj  is  sufficiently  small  we  may  expect  that  m1  =  m°  +  5m  will  more  nearly  satisfy  eq  56 
than  m°  did.  As  a  measure  of  a  model’s  suitability  ,  we  may  define 


5m.  =  D°  -  DAm°)  i  =  1,  . . .  K.  (59) 

,o  1  1  1 


Sm.  +0(|5m|2)  i  =  1,  . . .  K.  (58) 

,o  J 


((m)  =  £ 


\D.(m)  -  D°| 


i  =  l 


D; 


(60) 


Equations  59  constitute  a  K  x  3N  set  of  linear  equations  in  the  components  of  5m.  We  may 
not,  in  general,  expect  to  find  5m  exactly  satisfying  eq  59  for  all  possible  cases.  If  the  rank  of  the 
system  does  not  exceed  3 N,  such  a  5m  exists  but  is  not  necessarily  unique.  If  the  r^nk  exceeds 
3JV,  it  does  not  exist. 

Various  methods  of  solution  have  been  applied  to  the  system  59  (Backus  and  Gilbert,  1967; 
Anderson  and  Smith,  1968;  Smith  and  Franklin,  1969;  Jordan  and  Franklin,  manuscript,  1971).  We 
adopt  here  a  general  technique  proposed  by  Franklin,  (unpublished  manuscript,  1969.) 

We  rewrite  eq  59  more  compactly  as 
j4(m°)  •  5m  =  R(m°) 

where  4(m°)  is  the  matrix  whose  elements  are 


(61) 
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(62) 


and  R(m°)  is  the  vector  of  data  residuals  whose  ith  component  is  -  DjCm0).  We  now  regard 
eq  61  as  a  linear  relation  between  three  stochastic  processes:  a  signal  process,  a  data  process 
and  a  noise  process.  8m  is  a  sample  of  the  signal  process,  and  R(m°)  is  the  sum  of  a  sample  of 
the  data  process  plus  a  sample  of  the  noise  process.  Each  process  is  taken  to  have  zero 
expectation. 

The  use  of  stochastic  techniques  to  solve  eq  61  is  based,  in  part,  upon  contemplation  of  some 
of  the  potential  sources  of  error  entering  into  the  relation.  The  measured  spectral  values  D?  are 
contaminated  by  measurement  error  and  possible  mode  misidentification.  The  physical  system  upon 
which  measurements  are  made  may  deviate  from  the  class  of  models  M  in  which  m  must  lie.  It  is 
possible  that  there  is  no  model  m  in  M  that  would  then  satisfy  eq  56. 

Let  Rm  denote  the  autocorrelation  associated  with  the  signal  process  5m  ,  and  R°  denote  the 
autocorrelation  operator  associated  with  the  noise  process.  Rm  is  a  3 N  x  3 N  square  matrix  whose 
(i,  j)th  component  is  the  expectation  of  the  product  R®  is  defined  analogously.  The  best 

linear  estimate  for  Sm  is  t,iven  by  (Franklin,  1969): 


Sm  =  R°  •  i(m°)  •  U(m°)  •  Rm  •  4T(m°)  +  R0]-1  •  R(«°).  (63) 

A  solution,  5m,  can  be  guaranteed  to  exist  if  R°  is  a  positive  definite  matrix. 

Fending  the  acquisition  of  experience,  we  will  forego  at  this  time  any  suggestions  about  the 
fabrication  of  Rm  and  R°. 

We  have  outlined  above  a  method  for  inverting  measured  eigenvalue  spectra  to  produce  a  model 
m  consistent  with  the  spectra.  The  only  remaining  theoretical  problem  is  the  computation  of  the 
partial  derivatives, 


which  form  the  components  of  4.(m°).  These  expressions  are,  typically,  extracted  by  applying 
Rayleigh’s  principle.  Particularly  good  discussions  are  available  in  Backus  and  Gilbert  (1967) 
and  Dahlen  (1968).  We  give  here  only  the  results.  If  A,  p  and  p  in  a  particular  layer  are  altered  by 
small  amounts,  the  resulting  alteration  in  <u  is  given  by 

rN 

1/  (SAA  +  8pM  8pR)rzdr 1 

So)  -  -  (64) 

rN 

I2<u  /  [l/2  +  1(1  +  l)V2]przdrl 
0 

where 

(<?  rU  +  rlF )2, 


\ 


(65a) 
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M  =  2(<9r{/)2  +  1(1  +  1  )r2[U  +  (rdt-  1)F]2  + 

+  r_2F2  +  r2u  -  1)1(1  +  1  )(I  +  2 )U2, 

(65b) 

R  =  -ioZ[U2  +  HI  +  1)F2], 

(65c) 

F  =  2U  -  1(1  +  1)F. 

(65d) 

Preliminary  Experimental  Results 

As  a  first  step  in  the  validation  of  the  foregoing  theoretical  results  we  elected  to  attempt  to 
measure  the  normal  mode  spectra  of  layered  spheres  of  known  composition.  Spherical  samples  are, 
in  general,  difficult  to  prepare,  and  our  first  specimen  was  composed  of  a  thin  steel  shell  filled 
with  distilled  water. 

The  steel  shell  was  formed  by  joining  two  stainless  steel  hemispheres  with  epoxy.  The 
sphericity  and  homogeneity  of  the  resulting  shell  is  open  to  question;  the  hemispheres  were  pro¬ 
cured  from  an  industrial  float  works.  However,  we  believed  that  mechanical  perfection  was  not,  at 
this  point,  necessary. 

A  small  transducer  was  bonded  to  the  shell  and  the  sample  was  placed  in  the  experimental 
arrangement  shown  in  Figure  12.  Coupling  between  the  sample  and  the  "driving”  transducer  was 
weakened  by  inserting  rubber  padding  and  an  "0”  ring  between  the  two.  Some  experimentation 
indicated  that  coupling  in  this  configuration  is  sufficiently  weak  that  the  measured  spectra  are  not 
significantly  influenced. 

The  apparatus  produced  a  graphical  record  of  the  response  of  the  sample  to  a  sinusoidal  source 
as  a  function  of  frequency.  A  portion  of  that  result  is  recorded  in  Figure  13.  The  center  frequency 
of  significant  peaks  is  given  above  each  such  peak.  The  amplitude  scale  is  arbitrary,  but  linear. 
Tentative  assignments  of  peaks  to  a  particular  mode  are  indicated  by  the  designation  nSj,  where  n 
is  the  overtone  number  and  i  is  the  harmonic  degree. 

Figure  14  is  s.  juxtaposition  of  the  observed  resonance  frequencies  and  the  values  computed 
for  a  sphere  composed  of  an  inner  core  of  density  0.998,  having  a  eompressional  velocity  of 
1.498  x  10s,  and  a  radius  of  7.62.  The  core  is  surrounded  by  a  shell  of  thickness  0.0904,  with  a 
density  of  7.3772  and  having  eompressional  and  shear  velocities  of,  respectively,  5.79  x  10s  and 
3.1  10s.  We  regard  the  agreement  as  being  generally  satisfactory. 

In  Figure  13,  it  can  be  seen  that  we  have,  in  several  instances,  assigned  a  mode  to  a  small 
collection  of  adjacent  peaks.  Perturbation  theory  (Dahlen,  1968)  tells  us  that  small  deviations  from 
spherical  symmetry  will  generally  split  the  21  +  1  individual  modes  associated  with  a  given  (n,  1) 
pair  apart  in  frequency.  That  is,  asphericities  either  remove,  or  at  least  decrease,  the  degeneracy 
of  a  particular  mode.  (Only  the  family  nS0  is  not  degenerate,  and  it,  too,  can  still  be  shifted.)  It  is 
not  umeasonabJe  to  expect  that  the  abundant  asphericities  present  in  our  sample  will  produce  such 
an  effect  and  that,  for  instance,  the  twin  peaks  associated  with  0S2  in  Figure  13  are  an  expression 
of  this.  We  believe  that  such  effects  can  be  greatly  reduced  by  increasing  the  mechanical  symmetry 
of  the  sample. 

The  sample  fabrication  method  used  above  is  not  a  particularly  happy  one  for  routine  measure¬ 
ments  of  the  properties  of  fluids.  Figure  15  shows  the  partial  derivative  of  eigenfrequency  with 


Amplitude  (arbitrary  units) 
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Silicon*  Rubber 


Figure  12.  Experimental  apparatus  used  to  measure  frequencies  of  free  oscillation. 


Figure  13.  Portion  of  the  measured  spectrum  of  a  stainless  steel  shell  0.040  in.  thick  and 
6  in.  O.D.,  filled  with  distilled  water.  Peak  frequencies  are  shown  and  tentative  mode  assign¬ 
ments,  based  on  handbook  data,  are  indicated. 
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Figure  14.  Open  circles  represent  loci,  in  the  (I,  1) 
plane,  of  the  frequencies  of  free  oscillation,  computed 
for  a  stainless  steel  shell,  filled  with  distilled  water, 
of  tne  dimensions  used  in  the  experiment.  Filled  cir¬ 
cles  indicate  measured  resonances.  Arrows  indicate 
overlapping  symbols.  Assignment  of  an  ordinate  to 
measured  frequencies  is  arbitrary. 

respect  to  compressional  velocity  in  the  fluid  and  the  shell,  and  shear  velocity  in  the  shell,  as  a 
function  of  i,  for  the  fundamental  modes  oSj.  We  see  that  for  1  >  2  shear  velocity  in  the  shell 
dominates  all  other  controlling  parameters,  and  that  for  /  >  5  compressional  velocity  in  the  fluid 
is  the  least  significant  parameter.  At  1  =  13,  compressional  velocity  in  the  fluid  is  about  1/40 
as  important  as  shear  velocity  in  the  shell. 

A  preferable  arrangement  is  one  in  which  the  quantity  of  interest  is  the  controlling  parameter. 
One  way  to  achieve  this  is  to  utilize  higher-order  overtones.  For  example,  for  ^g  we  have 


o  Computed 
h  •  Observed 


°S 


•O 


—  =  0.126 


df 


1.65  v 


df 


3.39  x 


(fluid  core), 

10"^  (shell),  and 

10~3  (shell). 


Thus  jS3  is  substantially  more  influenced  by  the  properties  of  the  fluid  than  by  those  of  the  shell. 
This  behavior  results,  in  a  general  way,  from  the  increasing  concentration  of  energy  in  the  interior 
associated  with  higher  and  higher  overtones.  Unfortunately  one  result  of  this  is  that  such  modes 
are  difficult  to  excite,  or  observe,  from  the  surface. 
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Figure  IS.  The  partial  derivatives  of  the  fundamental  frequencies  of 
free  oscillation,  as  a  function  of  1  and  with  respect  to  various  indi¬ 
cated  properties,  of  a  water-filled  steel  shell.  For  visual  ease,  the 
derivatives  have  been  Joined  by  smooth  carves. 


A  more  promising  approach  is  to  use  shells  composed  of  a  material,  such  as  Lucite,  whose 
properties  do  not  represent  as  severe  an  impedance  contrast  to  the  interior  material.  Exploratory 
calculations  are  presently  being  made  on  the  optimum  shell  composition.  We  are  also  arranging  to 
have  shells  machined,  as  opposed  to  other  methods  of  fabrication,  to  greatly  reduce  the  variations 
and  asphericities  of  the  sample. 

The  preparation  of  spherical  soil  samples  proved  to  be  considerably  more  difficult,  owing  to  the 
volume  change  of  water  upon  freezing.  After  extensive  experimentation  we  were  finally  able  to  form 
a  frozen  soil  inside  a  spherical  steel  shell  by  maintaining  strong  temperature  gradients  across  the 
sample  while  allowing  for  the  venting  of  unfrozen  water  from  the  sphere. 

Figure  16  shows  a  portion  of  the  resonance  spectrum  for  a  4-in.-diameter  sphere  of  frozen,  fully 
saturated  20/30  Ottawa  banding  sand.  The  experimental  arrangement  used  in  this  measurement 
differed  from  that  depicted  in  Figure  12  in  that  a)  the  sphere  was  supported  on  inflated  plastic 
cushions  and  b)  the  applied  signal  consisted  of  a  four  cycle  tone  burst  generated  every  10  milli¬ 
seconds. 

The  modes  QSg  and  oS0,  together  with  their  partial  derivatives,  were  used  to  estimate  the 
compressional  and  shear  velocities  of  the  frozen  soil.  Forward  calculations  for  the  new  model  are 
shown  in  Table  IV  and  compared  to  the  data.  The  maximum  relative  error  for  the  lowest  five  ob¬ 
served  modes  is  3.7%  for  0Sg,  which  Figure  16  shows  as  being  badly  split 
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f  (kHz) 

Figure  16.  Portion  of  the  measured  spectrum  of  a  stainless  steel  shell,  0.040  in.  thick 
and  4  in.  O.D.  filled  with  frozen,  saturated,  20/30  Ottawa  banding  sand  at  -38nC. 

Table  IV.  Computed  and  observed  frequencies  of  tree  oscillation  of  a  saturated, 
frozen  20/30  Ottawa  banding  sand  sample  encased  in  a  steel  shell. 


Mode 

Computed' 

(kHz) 

Observed 

(kHz) 

Comp-Oba 

Comp 

oS, 

21.01 

21.01 

0 

oS, 

27.72 

28.11 

-1.4  x  102 

oS, 

30.96 

31.96 

-3.2  x  10* 2 

oSo 

33.42 

33.44 

-6  x  10*" 

oS4 

39.38 

39.88 

-1.3  x  10'* 

oSs 

47. 14 

?2 

1.  Computed  for  r,  =  4.961  cm,  V  .  =  4.391  x  10^  cm/sec, 

^1  = 

2.66  x  10°  cm/sec, 

r. 

Pj  =  2.0,  r2  =  5.062  cm,  K  g  = 

5.79  x 

10  cm/sec,  Kg«,  = 

3.1  x  10  cm/sec, 

p2  =  7.8772. 

2.  Splitting  for  oSg  is  too  severe  to  allow  a  useful  result. 

ELASTIC  PROPERTIES  OF  MATERIALS 


55 


Conclusions  and  Projected  Research 

We  believe  that  this  procedure  will  provide  a  theoretically  straightforward  and  experimentally 
practical  method  of  measuring  the  elastic  properties  of  solids  and  fluids,  and  in  particular,  of 
frozen  ground.  We  further  expect,  with  the  existing  apparatus,  to  be  able  to  measure  Q,  or  attenua¬ 
tion,  simultaneously  with  the  elastic  properties.  (The  necessary  theory  for  attenuation  measure¬ 
ments  is  well-developed.  See  Anderson  and  Archambeau  (1964).) 

The  technique  should  be  particularly  useful  in  discerning  small  variations  of  elastic  properties 
as  a  function  of  temperature.  Such  information  is  of  value  in  deciphering  the  physical  chemistry 
of  frozen  soils  through,  and  near,  the  freezing  point. 
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IV.  DETERMINATION  OF  ACOUSTIC  PROPERTIES  OF  FROZEN  EAKTH  MATERIALS 
BY  THE  USE  OF  A  CRITICAL  ANGLE  TANK 

by 

Y.  Nakano,  M.  Smith  and  R.  Martin 


Introduction 

One  of  the  standard  methods  for  determining  acoustic  properties  of  earth  materials  is  an 
ultrasonic  technique  using  a  critical  angle  tank.  In  spite  of  its  inherent  interpretive  problems,  it 
theoretically  allows  simultaneous  operation  in  both  dilatation  and  shear  at  wavelengths  small  com¬ 
pared  to  the  sample  thickness. 

The  method,  applied  to  metal,  was  first  described  by  Schneider  and  Burton  (1949)  and  sub¬ 
sequently  used  by  Subbarao  and  Rao  (1957)  on  rocks.  King  and  Fatt  (1962)  used  it  to  determine 
velocities  in  rock  specimens  subjected  to  confining  pressures,  and  Wyllie  et  al.  (1962),  Gregory 
(1963),  and  Banthia  et  al.  (1965)  applied  the  techniques  to  determine  velocities  in  rocks  under  both 
differential  external  confining  pressures  and  internal  pore  fluid  pressures.  Rec'mtly  Attewell  (1970) 
used  this  method  to  study  triaxial  anisotropy  of  rocks. 

In  the  present  work  dilatational  and  shoar  wave  velocities,  as  well  as  attenuation  of  several 
standard  frozen  soils,  were  determined  by  the  use  of  the  critical  angle  tank. 


Experimental  Apparatus 

The  critical  angle  tank  is  constructed  from  Vi- inch-thick  Plexiglas,  13  inches  long  x  9  inches 
wide  x  9  inches  deep  (Fig.  17).  At  either  end,  a  transducer  mounted  on  an  aluminum  rod  is  inserted 
through  the  wall  of  the  tank  by  means  of  an  O-ring  seal.  The  distance  between  the  two  transducers 
can  be  changed.  One  transducer,  made  of  2.5-inch-diaraeter  x  0.5-inch-thick  PZT4,  serves  as  a 
transmitter  for  producing  uniform  plane  dilatational  waves.  The  transmitter  is  connected  to  a  Wayne 
Xerr  Model  SR-268  signal  generator  via  a  General  Radio  Model  1396-B  tone-burst  generator  so  that 
a  harmonic  burst  of  any  number  of  cycles  can  be  sent  to  the  transmitter  with  any  desired  time  interval. 
The  other  transducer,  of  1.0-inch-  diameter  x  0.5-inch-thick  PZT4,  serves  as  a  receiver,  and  is 
connected  to  a  Tektronix  type  567  dual-trace  oscilloscope  via  a  Krohn-hite  filter,  which  passes  only 
signals  of  a  specified  frequency  range.  One  trace  of  the  oscilloscope  is  triggered  directly  from  the 
pulse  generator  to  display  received  signals.  The  other  oscilloscope  trace  is  used  for  a  Rutherford 
Electronics  time  delay  generator,  which  is  triggered  by  the  pulse  generator  and  is  used  for  accurate 
measurements  of  arrival  time  with  error  less  than  +0.1  sec. 

The  tank  is  filled  with  Dow  Corning  200  silicone  oil  used  both  as  a  cooling  medium  and  as  a 
medium  for  transmitting  the  input  pulses  across  the  tank  to  the  receiver.  The  cooling  medium  is 
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Figure  17.  Critical  angle  tank. 


circulated  through  a  Forma  Scientific  Model  2095  refrigerated  bath  which  maintains  the  required 
temperature.  A  sample  holder  7  inches  long  x  8  inches  deep  that  can  hold  a  specimen  up  to  1  inch 
thick  is  placed  vertically  inside  the  tank.  A  frozen  soil  specimen,  4.5  inches  in  diameter,  is  mounted 
at  the  center  of  the  sample  holder  so  as  to  intercept  the  ultrasonic  beam.  The  surface  of  the  frozen 
soil  specimen  is  sealed  by  a  thin  vinyl  sheet.  The  sample  can  be  rotated  about  a  vertical  axis 
running  through  the  center  of  the  tank,  the  angle  of  rotation  of  the  sample  from  normal  beam  incidence 
being  indicated  on  a  circular  dial  graduated  in  degrees  and  readable  to  0.1°  by  means  of  a  vernier 
scale. 


Theory 


A  transmitted  sound  wave,  after  passing  through  the  fluid,  strikes  the  sample  at  an  angle  i  to 
its  normal  (Fig.  18a),  Of  the  incident  dllatational  energy,  part  is  reflected  back  at  an  angle  i  and 
part  is  transmitted  into  the  sample.  Since  the  dilatational  velocity  Cp  in  frozen  soils  (~  4.0  km/sec) 
greatly  exceeds  that  (Cw)  in  silicone  oil  (~  1.0  km/sec),  the  transmitted  (refracted)  portion  of  the 
incident  energy  is  rotated  from  the  normal  at  an  angle  0p  >  i  according  to  Snell’s  law: 


(1) 
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Figure  18.  Reflection  and  refraction  geometry  of  a  beam  of  dilatational  pulses 
on  a  slab  of  frozen  soils  in  silicone  oil.  (a)  Near  normal  incidence,  (b)  Crit¬ 
ical  refraction  of  the  dilatational  beam,  (c)  Critical  refraction  of  the  shear 

beam. 


The  component  of  the  particle  motion  not  accounted  for  by  the  transmitted  dilatation  P  forms 
a  shear  wave  S  which  is  refracted  at  a  small  angle  0g  ( <0p )  since  Cp  >  Cg  (~2.0  km/sec)  >  Cw. 
S  is  polarized  within  the  plane  of  incidence  again  according  to  Snell's  law: 


At  the  second  surface  of  the  sample,  since  fluid  cannot  accommodate  the  shear  mode,  both  wave- 
trains  are  refracted  back  towards  the  normal  at  an  angle  i  in  the  form  of  dilatational  waves  and 
activate  the  receiving  transducer.  Each  of  the  shear  and  dilatational  signals  striking  the  second 
surface  will,  in  general,  give  rise  to  both  shear  and  dilatational  waves,  in  the  sample,  which 
emanate  from  the  interface.  However,  by  considering  only  the  first  portion  of  the  received  signal 
we  may  eliminate  these. 

As  i  is  increased  by  mechanical  rotation  of  the  sample,  a  critical  angle  i  is  reached  (Fig. 
18b)  at  which  P  skims  across  the  surface  of  the  sample  in  a  transitional  stage  between  refraction 
and  reflection.  At  this  stage,  from  eq  1, 


Cw 

8<n  (icp)  ‘ 


(3) 


As  i  is  increased  still  farther,  P  becomes  totally  reflected  and  at  a  second  critical  angle,  icg 
(>  icp),  S  skims  across  the  surface  of  the  sample  (Fig.  18c).  (In  fact,  reflection  is  complete  only 
if  the  sample  is  infinite  in  thickness.  Signals  incident  at,  or  beyond,  thy  critical  ar "le  will  ‘‘tunnel" 
through  a  finite  sample  in  a  fashion  exactly  analogous  to  the  well-know  i  quantum  mechanical 
phenomenon.  For  samples  greater  than  a  few  wavelengths,  we  may  disrtgard  this  effect.)  Again, 
from  eq  2, 


C„  = 


sin  0„> 


(4) 
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At  normal  incidence  the  beam  travels  the  shortest  dimension  of  the  sample  and  suffers  minimal 
attenuation.  As  i  is  increased,  the  wave  path  in  the  sample  increases  and  the  transmitted  pulse 
amplitude  decreases.  When  i  =  icp,  P  is  removed  from  the  transmitted  beam  and  the  resolved 
amplitude  for  i  exhibits,  in  practice,  an  intermediate  minimum  rather  than  an  inflection  on  the 
curve.  At  i  =  i  ,  S  is  eliminated  from  transmission  and  one  is  left  with  an  absolute  minimum 
amplitude  In  practice,  the  theoretical  sequence  of  events  outlined  above  is  disorganized  as  a  re¬ 
sult  of  the  interference  between  the  incident  beam  and  reflected  beams  within  the  specimen.  Inter¬ 
ference  patterns  often  lead  to  difficulty  in  locating  icp  with  precision.  Therefore,  another  method, 
a  fluid  displacement  technique,  was  used  to  determine  dilatational  velocities  in  the  present  work. 
The  difference  in  delay  time  St  between  the  pulses  received  after  crossing  the  liquid  path  alone 
as  compared  with  their  passage  through  liquid  and  specimen  interposed  at  normal  incidence  produces 
a  value  for  C  in  the  following  way.  If  tj  ia  tuo  time  delay  arising  from  transmission  across  the 
tank  through  liquid  alone  and  iw  is  the  wave  path  in  liquid,  then  tj  =  iw/Cw.  The  time  delay 
arising  from  transmission  across  the  tank  when  a  sample  is  interposed  at  normal  incidence  is, 


where  d  is  the  sample  thickness.  Then  the  difference  in  both  delays  is, 


The  attenuation  is  determined  from  the  measurement  of  the  ratio  of  the  amplitudes  Aj  and  A. 
of  the  received  waves  for  specimens  of  two  different  thicknesses  Li  and  L-,  respectively 
(Auberger  and  Rinehart,  1961).  Namely,  the  attenuation  coefficient  a  (nepers/ cm)  is  given  as: 


In  (Aj/Aj) 


(Li  -  Li> 

or  the  value  of  Q  is  given  as: 

Or/ A)  (Ls  -  LA 

Q  =  - - - - 

In  (A,/Aj) 


where  A  is  a  wavelength. 
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Results  sad  Discussion 

We  examined  three  standard  types  of  soils,  Ottawa  sand,  Hanover  silt,  and  Goodrich  clay,  whose 
gradation  curves  obtained  according  to  ASTM  test  procedures  are  shown  in  Pig.  1.  The  size  of  the 
specimen  used  for  velocity  measurements  is  4.5  inches  in  diameter  x  1.0  inch  thick.  All  measure¬ 
ments  were  made  at  subzero  temperature  due  to  the  fragility  of  the  sample  holding  device. 

Dilatational  wave  velocity 

Velocity  of  dilatational  waves  was  measured  as  a  function  of  temperature  and  frequency.  It 
was  found  that  the  frequency  has  no  distinguishable  effect  on  dilatational  velocities  in  the  range 
between  300  kHz  and  1.2  MHz  examined.  Fig.  19  shows  the  observed  dilatational  velocities  as  a 
function  of  temperature.  The  detailed  discussion  on  the  correlation  between  the  dilatational 
velocity  and  temperature  or  unfrozen  water  content  is  presented  in  another  part  of  this  report. 


Tamptraturt,  C 

Figure  19.  Dilatational  velocity  vs  temperature. 


Shear  wave  velocity 

Measured  shear  velocities  are  shown  in  Figure  20  as  a  function  of  temperature.  Although  the 
data  are  somewhat  scattered,  there  is  a  general  tendency  for  the  shear  velocity  to  decrease  with 
ascending  temperature.  The  roll  of  unfrozen  water  in  shear  wave  propagation  is  more  complicated 
than  the  decrease  of  velocity  as  in  dilatational  wave  propagation.  Since  liquid  water  cannot 
accommodate  the  shear  mode,  shear  waves  are  expected  to  attenuate  severely  while  traveling  through 
a  layer  of  liquid  water.  Therefore  the  crystalline  matrices  consisting  of  soil  minerals  and  possibly 
ice  are  the  path  through  which  shear  waves  propagate.  In  view  of  the  fact  that  the  slope  of  velocity- 
temperature  curves  does  not  differ  markedly  among  the  tested  samples,  soil  mineral  matrix  might 
play  a  major  roll.  Since  the  soil  samples  used  are  well  packed,  soil  minerals  must  contact  each  other 
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Figure  20.  Shear  velocity  vs  temperature  for  20-30  Ot¬ 
tawa  sand,  Hanover  silt,  and  Goodrich  clay. 


closely.  If  the  soil  is  not  well  packed,  a  different  result  might  be  obtained.  The  shear  velocities  of 
crystalline  rock  and  polycrystalline  ice  are  about  3.0  km/sec  and  1.6  km/ sec  respectively.  The 
measured  shear  velocities  of  frozen  soils  fall  between  these  two  bounds. 

The  shear  velocity  was  found  sensitive  to  frequency  in  an  almost  linear  manner  in  the  range  of 
0.3  MHz  to  1.2  MHz.  The  velocity  increases  with  increasing  frequency.  This  trend  was  also  re¬ 
ported  by  Attewell  (1970),  who  measured  the  shear  velocity  of  the  Penrhyn  slate  of  North  Wales  using 
the  critical  angle  method.  He  claimed  that  the  shear  wave  dispersion  resulting  in  frequency  de¬ 
pendence  of  velocity  is  a  function  of  the  critical  angle  technique  and  does  not  represent  the  shear 
velocity  conditions  in  the  slate.  Any  increase  of  velocity  with  frequency  implies  distortion  of  the 
spectrum  of  traversing  waves  through  the  solid.  Although  Attewell’s  reasoning  based  upon  much 
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experimental  evidence  is  convincing,  the  reason  why  the  critical  angle  technique  does  not  accurately 
represent  the  shear  conditions  in  the  slate  was  not  explained. ' 

Since  we  do  not  have  sufficient  experimental  results,  we  are  not  able  to  evaluate  the  accuracy 
of  shear  velocity  measurements  in  frozen  soils.  Further  efforts  are  needed  to  clarify  this  problem. 

Attenuation  of  dilatational  waves 

Attenuation  of  waves  in  earth  materials  at  high  frequency  is  one  of  the  least  understood 
phenomena.  Wyllie  et  al.  (1963)  stated  that  attempts  to  measure  attenuation  in  porous  media  at 
frequencies  higher  than  100  kHz  gave  equivocal  results.  Auberger  and  Rinehart  (1961)  measured 
attenuation  of  longitudinal  waves  in  rocks  by  the  pulse  method  in  the  frequency  range  0.3  MHz  to 
1.0  MHz.  No  definite  conclusion  was  made  on  the  correlation  between  frequency  and  attenuation, 
since  attenuation  does  not  follow  any  marked  law  of  increase  or  decrease  with  frequency.  Attewell 
(1970)  reported  attenuation  of  dilatational  waves  in  hard  blue  Penrhyn  slate  by  the  critical  angle 
technique  in  the  frequency  range  0.5  MHz  to  5.0  MHz.  He  did  not  describe  any  difficulty  or  problem 
concerning  attenuation  measurement. 

We  measured  the  attenuation  of  dilatational  waves  in  frozen  soils  by  using  samples  of  different 
thickness.  Although  the  wave  patterns  of  the  received  signals  for  samples  of  different  thickness 
resemble  each  other  and  the  amplitude  of  an  initial  rise  or  first  arrival  can  be  measured  without 
difficulty,  we  are  plagued  by  lack  of  reproducibility.  We  do  not  know  whether  this  is  inherent  in 
the  frozen  soils  examined  or  is  due  to  error  in  measurements.  Therefore  we  are  not  able  to  present 
our  results  in  this  report. 
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APPENDIX  A.  ANALYTIC  SOLUTIONS  OF  THE  TRANSFORMED  RAVE  EQUATION 

We  represent  u  by 

u  =  V^  +  y  x  Try  +  V  x  V  x  Tra  (Al) 

where  <p ,  y  and  a  are  scalar  fields.  The  representation  Al  was  chosen  because  it  is  "natural”  to 
the  field  equations  (eq  3). 

It  is  helpful  to  develop  the  following  useful  expressions 

V  •  u  =  V2^', 
and 

VxVxu  =  yxVx  ^  x  (Try)  +  V  x  V  x  V  x  Vx  (Fra)  = 


=  V  x  lv(V  •  Try)  -  V2(rrx) I  + 

+  V  x  V  x  lv(V  •  Tra)  -  V2(rr<7)i . 

We  expand  VzYry)  as 

V*Vrx)  =r|v2(rx)  -  +  Vjj^xj 

=  Tr  v2^  +  Vl2x  ! 

as  may  easily  be  shown  by  expanding.  So, 

yxyxu  =  -yxTr  v2X  -  V  x  V  x  Tr  V2<r. 

If  we  insert  these  into  the  field  equation  (eq  3)  and  regroup  terms,  we  find 

Vi pd^xjj  -  (A  +  2p)V2^!  +  V  *Tr  \pdZy  -  /iy2^i  + 

+  Vx  Vxrripd2c7-/iy2c7|  =  0. 

In  order  to  ensure  that  eq  3  is  satisfied,  it  is  sufficient  that  ,  y  and  a  be  solutions  to 
pdzif/  -  (A  +  2  p)v2|A  =  0, 


(A2a) 


PtfKDMG  PAGE  BUNK 
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P^X  ~  PVZX  =  0,  (A2b) 

and 

pdfo  -  pVZa  =  0.  (A2c) 

We  have  not  shown  that  all  solutions  to  eq  3  can  be  expressed  in  terms  of  functions  satisfying 
eq  A2  and  Al.  These  solutions  are,  however,  known  to  be  complete  (Sternberg,  1960). 

We  now  Fourier  transform  the  system  and  introduce  the. expansions 

*(r,  0,</>,<o)  =1  {  0f  (r,  to)  Ff  (0,  <f>),  (A3a) 

1=0  m=-l 

x(r,  0,  &co)  =  l  I  yf  (r,  <o)Y?  (0.  0),  (A3b) 

1=1  m=-l 

and 

air.  6,  <f>,  co)  =  1  £  of  (r,  <u)  Ff  (0,  $).  (A3c) 

1=1  m=-l 

The  terms  of  degree  1  =  0  in  the  expansions  for  y  and  a  have  been  omitted  since  they  do  not 
contribute  to  the  displacement  field. 

The  expansions  A3a-A3c  are  inserted  into  the  transformed  versions  of  A2a-A2c.  We  make  use 
of  eq  6  and  14  to  simplify  the  result.  If  /  is  some  scalar  field,  then  by  eq  6 

VI  =Tdtf  +  Vi('~M 


and 


V2/  =  d*l  +  |  0tl  +  i  i 


-2_2 


(A4) 


The  resulting  expressions  are  multiplied  by  F™  (0,  <f>)  sin  6  and  integrated  over  6  and  <£ .  We  appeal 
to  the  orthogonality  relation  12.  We  then  have 


L2  +  !  _io_Li)i  r 

|  r  r  r  A  +  2fi  f2  p  ’ 

J,2  2  ,  <jzp  1(1  +  1)1  m, 

V-  *  *  7 - — /  *■  <'• 


A  =  0, 


)  =  o. 


(A5a) 


(A5b) 


and 
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i(i 


^}«r  <'• 


<u) 


0. 


(A5c) 


Each  of  the  operators  in  brackets  is  some  form  of  the  spherical  Bessel’s  operator.  Its  solutions 
are  well-known  and  they  are 


and 


where 


and 


■T 


k  = 


=  -4™  (o>) /j  (Irr)  +  Bf  fo>) yj  Ur), 

(A6a) 

=  Cf  (o>)  /j  (yr)  +  D™  (<u)  yj  (yr). 

(A6b) 

=  Ef  (oJijtyr)  +  Ff  (a>)y,(yr). 

(A6c) 

6) 

/A  +  2/i 

(A7) 

V  ' 

Cl) 

—  • 

(A8) 

OJ 

Y  =  -r=-  • 

jf 


We  wish  now  to  relate  U,  V  and  V/  to  \  and  a.  To  do  this,  we  will  rearrange  eq  A1  to 
resemble  eq  4.  For  convenience  we  will  drop  subscripts  and  superscripts. 

We  note  immediately  that 


V0  =  Tdt<jj  + 


Also, 


V  x  ~&X  =  rXV  *  r-Tx  v(ry) 
=  -T x  V(rx) 


•  -  r  x  VjX. 

The  third  term  can  be  expanded  as 

V  x  v  x  Tro  =  V(V  •  Tro)  -  Vs  {rro) 


=  -7rv8o  -  v(2o)  +  V  [r_8<?r  (r3®)] 
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V  x  V  x  Tta  =  Tl-rv2^  +  d,  [t~Zdth3o)  -  Za]\  + 
+  Vi  ir“3dI(rso)  -  2t~la] 

=  Tir1v2ffl  +  Vi  ir_1a  +  dra I 


=  7^-  ^  ^  +  Vt  ir_1<?r  (ro)  I. 

Collecting  these  we  have 

u  =T^dt<//  -  — 1  V.  a ^  +  Vt  ir-10  +  r_1(9r(ra)|  -  r  x  V^. 
Therefore 

u®  =  dr4>y  -  — -t-  -  o“, 

Vf  =  r1^  +  ar(rof)]. 


(A9a) 

(A9b) 
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APPENDIX  B.  THE  TRACTION  ON  A  SPHERICAL  SURFACE 

Let  X  be  the  elastic  stress  tensor  resulting  from  a  displacement  field  u.  T  is  given  by 

T  =  A(V  •  u)J  +  n(Vn  +  uv)  (Bl) 

where  uV  is  simply  the  transpose  of  Vu.  The  force  (not  stress)  acting  across  a  surface  whose 
unit  normal  islfis  given  by  TP  •  T,  which  is  a  vector  quantity.  We  represent  the  force  acting 
across  a  surface  whose  normal  is  the  radius  vectorTby 

T •  T  =TP  +  Vi<?  -  Tx  Vj R.  (B2) 

Our  problem  is  to  relate  the  scalars  P,  Q  and  R  to  the  scalars  U,  V  and  W  which  characterize  the 
displacement  field. 

We  note  first  that 

r* •  T  =  A(V  •  u)T  +  /rr* •  f  Vu  +  uVi .  (B3) 

The  divergence  of  u  can  be  expanded  as 

V  •  u  =  (T<5r  +  r~lVi)  •  Vu  +  VjF)  (B4) 

since 

V-  (-Tx  HO  -  -F.  (V  x  VjHO 
=  -7* •  (Tx  V<?rW) 

=  0 

because?*  x  VdfW  must  be  normal  toT.  Equation  B4  can  be  written  as 

V  .  u  =  dfU  +  r-1V \V  +7*^  •  VjF  +  r-1  Vj  •  71/ .  (B5) 

The  third  term  vanishes  since  df  commutes  with  V1(  and  Vi<?rv  is  normal  to 7.  The  fourth  term  is 
equal  to  (2/r)  U  as  may  be  seen  by  replacing  r-1  Vt  with  V  -  7<9r,  an  equivalent  expression.  We 
then  have 

V  .  n  =  (ar  +  j)u  +  r_l V? v«  (B6) 


,s 

$ 


70 


APPENDIX  B 


The  term  7*  •  I V  u  +  u  V  i  in  eq  B3  is  more  difficult  to  deal  with.  In  terms  of  the  coordinate 
representation  of  u,  (i.e.  uf,  ug,  and  u^), 

T»  iVu  +  u  Vi  =  7(2dfuf)  +  +  tdf 

'  ?(nsnr  V.  *  ■>'  (^)l  •  <B7» 

By  inspection  of  eq  5  which  explicitly  gives  the  coordinate  components  of  Vlt  we  see  that  eq  B7 
may  be  ro  writ  ten  as 


r  • 


I  Vu  +  u  Vi  =  r*(2 <?r  ur)  + 


+ 


We  note  that  dr  commutes  with  6  and  <f>.  The  expression  in  square  brackets  represents  the  non-radial 
portion  of  u  and  must  therefore  be  identical  with  Vt  v  -  1*  x  Vj  w.  Since  ur  is  identical  with  U, 
we  have 


or 


I  Vu  +  nVi  =  "r^2<5fv) 


I  Vu  +  u Vi  =  7t2dfu) 


VjU  +  r<Jr 


1-* 

—  r  x 

r 


Combining  the  above  results,  we  have 


(B8a) 

(B8b) 


(B8c) 


71 


APPENDIX  C.  NUMERICAL  TECHNIQUES 

We  outline,  briefly,  in  this  appendix  die  numerical  techniques  used  in  this  study  to  a)  generate 
the  suite  of  normal  modes  for  a  layered  elastic  sphere  and  b)  evaluate  various  integrals  of  interest 
associated  with  these  modes. 

For  a  given  1  and  a>,  the  generation  of  solution  functions  proceeds  exactly  as  outlined  in 
Section  B.  The  matrix  inversion  required  by  eq  45  and  46  was  not  done  explicitly.  We  chose  instead 
to  solve  two  sets  of  simultaneous  linear  equations.  The  Crout  Reduction  (Hildebrand,  1956)  was 
found  to  be  particularly  convenient. 

Bessel  functions  were  generated  by  using  Miller’s  well-known  recurrence  algorithm  (Abramowitz 
and  Stegun,  1968).  One  consequence  of  this  technique  is  that  the  accurate  evaluation  of  a  spherical 
Bessel  function  for  many  values  of  its  argument,  as  required,  say,  for  integration,  is  a  time-consuming 
process.  For  such  applications  it  would  porhaps  be  more  efficient  to  numerically  solve  Bessel’s 
equation,  but  computer  memory  limitations  did  not  permit  the  additional  coding  this  required. 

In  practice,  the  program  was  assigned  a  model  and  a  value  of  1  and  proceeded  to  compute  trail 
solutions  for  evenly  spaced  values  of  frequency.  As  the  computation  proceeded,  indicator  variables, 
as  explained  in  Section  B,  were  monitored  for  a  change  of  sign,  which  was  taken  to  indicate  a  zero 
crossing.  When  this  occurred,  an  estimate  was  made  of  the  location  of  the  zero  crossing  and  the 
algorithm  described  below  was  invoked  to  iteratively  improve  the  estimate.  In  general,  two  applica¬ 
tions  of  the  following  procedure  sufficed  to  locate  the  eigenfrequency  to  within  one  part  in  10*. 

Gilbert  and  Backus  (1967)  observed  that  Rayleigh’s  principle  could  be  utilized  to  improve 
estimates  of  eigenflrequencies  obtained  by  coarser  methods.  Suppose  that  for  some  frequency  «, 
near  an  eigenvalue,  cj*,  we  have  computed  a  trial  solution  and  find  that  the  stress-flree  surface  con¬ 
dition  cannot  be  met.  We  may  apply  Rayleigh's  principle,  or  perturbation  theory,  to  the  solution  we 
have  generated  to  estimate  the  change  in  <u  the  elimination  of  surface  stress  would  produce.  The 
first  order  estimate  for  this  change  is  given  by 

dw  =  (8r<u)  - .  (Cl) 

/  przWf(r)  +  1(1  +  l)Ff(r)ldr 

o 

We  will  not  derive  eq  Cl  here.  We  then  replace  w  by  u  +  5cj  and  repeat  the  process.  We  chose  to 
terminate  the  iteration  when  \5<o/(8co  +  <u)|  fell  below  1(T*. 

The  last  point  we  wish  to  mention  is  the  evaluation  of  integrals  of  the  form 

/  =  J  Z(n)dr  (C2) 

0 
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where  Z  is  some  operator  on  the  solution  functions,  U(r),  etc.  In  general,  we  may  expect  Z  to  vary 
appreciably  over  a  length  L  of  about 

2nVn 

L  =  - P  (C3) 

a 

where  Vp  is  the  local  congressional  velocity.  L  is  simply  the  wavelength  of  a  compressional  wave 
of  angular  frequency  a>.  In  computing  /,  the  program  was  designed  to  utilize  steps  not  exceeding 
<Lt;  where  «  is  a  small  (~  3  x  10*1)  number  and  Li  is  the  scale  length  appropriate  to  a  given  shell. 
This  technique  yielded  a  reliably  constant  accuracy  over  many  wide  variations  of  scale  without 
extracting  undue  computing  labor  for  small  values  of  <u. 


